Au niveau de la lisibilité, je trouve préférable d'écrire
momentum = 0.0d0
do i = 1, size(elements,1)
momentum = momentum + calculateMomentum(elements(i))
end do
plus de
momentum = sum( [( calculateMomentum(elements(i)), i=1, size(elements,1) )] )
car la première version a la forme de définir la valeur de l' momentum , tandis que la seconde correspond à une instruction d'accumulation de plus bas niveau. La différence devient plus prononcée dans le code du monde réel avec plus d'arguments et éventuellement plusieurs indices.
Cependant, la première version alloue un tableau temporaire. En tant que programmeur humain, je sais que cela pourrait être optimisé, alors je me demandais si Fortran propose une syntaxe, qui permet de calculer la somme sans boucle explicite ni tableau temporaire.
Pas de telle construction. Il semble qu'il n'y ait pas de syntaxe, comme ce que je demandais. Voir la réponse de Vladimir F.
C'est important, mais moins que je ne le pensais. J'ai fait mon propre benchmark ( pastebin , embarqué ) en utilisant la multiplication matricielle en utilisant plusieurs variantes.
C(i,j) = C(i,j) + A(i,k) * B(k,j) était le plus lent, probablement en raison de l'accès inutile au tableau à chaque étape.
C(i,j) = sum ([(A (i, k) * B (k, j), k = 1, N)]) `était environ 10-20% plus rapide que (1), malgré le tableau temporaire .
tmp = tmp + A(i,k) * B(k,j) , c'est-à-dire en utilisant une variable d'accumulation temporaire, était environ 20% plus rapide que (2).
C = matmul(A,B) était de loin le plus rapide, comparé à (3) par un facteur de 25 pour les matrices 500x500, passant à 50 pour 2000x2000, tandis que la vitesse relative des autres variantes est restée à peu près la même.
Conclusion: lorsque la tâche ne peut pas être exprimée facilement dans des fonctions de bibliothèque ou intrinsèques optimisées, la variante de sum a des performances viables et ne doit être optimisée que si les performances comptent vraiment à un tel degré dans cette partie du code.
3 Réponses :
Le fait qu'un tableau temporaire soit alloué ou net dépend des optimisations du compilateur. L'allocation de pile est de toute façon presque gratuite. La copie des valeurs prendra probablement plus de temps.
Le compilateur peut optimiser les étapes inutiles s'il peut s'assurer que le résultat sera le même. Cependant, il n'y a pas de syntaxe spéciale pour cela. Fortran essaie généralement de rester loin de l'implémentation réelle et laisse beaucoup sur le compilateur.
Il y a un constructeur de tableau dans l'instruction, mais cela peut être évité.
En Python, il y aurait sum([f(x) for x in xs]) vs sum(f(x) for x in xs) . La première version crée une liste inutile, comme Fortran crée un tableau inutile. J'espérais en quelque sorte m'assurer que l'intrinsèque puisse faire quelque chose de plus similaire à r = 0; do i=1,size(xs); r = r + f(xs(i)); end do tout en conservant le style de code plus déclaratif.
@kdb Il n'y a rien de tel à Fortran. Pas d'itérateurs, pas de fermetures, pas d'évaluation paresseuse.
Techniquement, Fortran ne crée pas de tableau inutile. C'est le compilateur que vous utilisez qui crée le tableau temporaire. De plus, au lieu de deviner ce que fait gfortran, vous pouvez regarder les décharges intermédiaires. Essayez -fdump-tree-original et -fdump-tree-optimisé (le d de fin peut être incorrect). Enfin, au lieu d'essayer d'être intelligent, écrivez le code de la manière que vous désirez, c'est-à-dire r = 0, do i=1, size(xs); r = r + f(xs(i)); end do
@evets Peut-être vouliez-vous adresser cela à kdb?
Pour expérimenter, j'ai essayé ce code (qui calcule la somme de l'inverse de arr ).
val = 0
do i = 1, num
val = val + testfunc( arr(i) )
enddo
Ensuite, sur mon ordinateur (mac2012), "gfortran-10 -O2 test.f90 && time ./a.out" donne
val = 18.997896413852555 in 0.225999996 (s) val = 18.997896413852555 in 0.252000004 (s) val = 18.997896413852555 in 0.246999994 (s) ... real 0m3.909s user 0m3.009s sys 0m0.867s
et "gfortran-10 -O3" donne
val = 18.997896413852555 in 0.437000006 (s) val = 18.997896413852555 in 0.453999996 (s) val = 18.997896413852555 in 0.437999994 (s) ... real 0m5.946s user 0m5.069s sys 0m0.842s
Dans les deux cas, la commande htop affiche ~ 1,5 Go alloué, ce qui peut être raisonnable si les lignes L1 et L2 utilisent une matrice temporaire (chacune ~ 800 Mo avec ~ 0,3 s pour l'allocation).
Comme il n'y a pas de syntaxe pour créer des "itérateurs", j'ai essayé de faire de testfunc() un élémentaire (ou un élémentaire impur). La seule différence ici est les lignes marquées avec <-- .
program main
use iso_fortran_env, only: dp => real64
implicit none
real(dp) val
real(dp), allocatable :: arr(:)
integer num, loop, i, t1, t2, trate
num = 10**8
arr = [( i, i = 1, num )]
do loop = 1, 10
call system_clock( t1 )
val = sum( testfunc( arr ) ) !<--
call system_clock( t2, trate )
print *, "val = ", val, " in ", (t2 - t1) / real(trate), " (s)"
enddo
contains
impure elemental & !<--
function testfunc( x ) result( ret )
real(dp), intent(in) :: x
real(dp) :: ret
ret = 1.0_dp / x
end
end program
Ensuite, "gfortran-10 -O2" donne
val = 18.997896413852555 in 0.875000000 (s) val = 18.997896413852555 in 0.888000011 (s) val = 18.997896413852555 in 0.833000004 (s) ... real 0m9.986s user 0m5.738s sys 0m4.210s
et "gfortran-10 -O3" donne
val = 18.997896413852555 in 1.02999997 (s) val = 18.997896413852555 in 1.10099995 (s) val = 18.997896413852555 in 1.17600000 (s) ... real 0m12.575s user 0m8.142s sys 0m4.387s
La commande htop affiche ~ 800 Mo, il semble donc que seul arr est alloué.
À titre de comparaison, le code suivant calcule val avec une boucle do explicite (en utilisant une version scalaire de testfunc() )
program main
use iso_fortran_env, only: dp => real64
implicit none
real(dp) val
real(dp), allocatable :: arr(:)
integer num, loop, i, t1, t2, trate
num = 10**8
arr = [( i, i = 1, num )] !! L1
do loop = 1, 10
call system_clock( t1 )
val = sum( [( testfunc( arr(i) ), i = 1, num )] ) !! L2
call system_clock( t2, trate )
print *, "val = ", val, " in ", (t2 - t1) / real(trate), " (s)"
enddo
contains
function testfunc( x ) result( ret )
real(dp), intent(in) :: x
real(dp) :: ret
ret = 1.0_dp / x
end
end program
qui a donné le même timing avec le deuxième code avec elemental + testfunc(arr) ci-dessus (sur mon mac).
Dans tout le code ci-dessus, l'option -Ofast entraîné une erreur d'exécution ( Illegal instruction ). Mais cela était dû à la ligne L1 ( arr = [( i, i = 1, num )] ). Si -Ofast arr avance et que je le remplis avec une boucle explicite, -Ofast également fonctionné sans problème (donnant presque le même timing avec -O3 dans ce cas).
Donc quelle est la réponse?
peut-être "il n'y a pas d'itérateur dans la norme" ...?
Alors pourquoi n'avez-vous pas mis ce i dans votre réponse? Peut-être que je suis juste en train de râler et que les gens votent pour votre réponse et ne votent pas pour la mienne, mais à mon avis, la réponse devrait d'abord et avant tout répondre à la question posée. Et puis montrez des expériences intéressantes sur les performances. Peut-être que je suis juste od et grincheux.
@VladimirF Salut, je pensais juste que ta (première) réponse donnait déjà la réponse rigoureuse (de la norme), ce qui, je pense, devrait être acceptée (cc: @kdd), j'ai donc ajouté quelques expériences sous un angle différent (car OP semblait s'intéresser aussi aux codes connexes). Quoi qu'il en soit, je pense que "élémentaire" est une fonctionnalité très utile.
OK, je pense qu'il est toujours bon de dire le point saillant d'une nouvelle réponse, même si cela a déjà été dit.
Si vous créez la fonction calculateMomentum et élémentaire, elle peut être utilisée à la fois pour les valeurs scalaires et pour les tableaux
par exemple avec:
type body
real :: mass
end type
vous pouvez appliquer ce qui précède à un tableau de v
integer::i type(body) :: ball real, allocatable :: v(:), p(:) real :: tot_p allocate(v(10)) v = [ (10+i, i=1, 10) ] p = calculateMomentum(ball, v) tot_p = sum(p)
avoir un tableau intermédiaire pour contenir les valeurs est avantageux car il maintient les données à proximité (probablement dans la ligne de cache) et la fonction sum() serait aussi rapide que possible.
imaginez le body du type comme suit par exemple
elemental function calculateMomentum(obj, v) result(p)
class(body), intent(in) :: obj
real, intent(in) :: v
real :: p
p = obj%mass * v
end function