Optimisation et Sommaire

parallélisme avec OpenMP FI 3/2007

l’approche grain fin 1 Optimisation et parallélisme
Martin.Jucker@epfl.ch, EPFL-SB-CRPP, avec OpenMP:
Vincent.Keller@epfl.ch, EPFL-STI-ISE-LIN, l'approche grain fin
Gonçalo.Pena@epfl.ch, EPFL-SB-IACS-CMCS et Martin Jucker,
Riccardo.Puragliesi@epfl.ch, Paul Scherrer Institut (PSI) Vincent Keller,
Gonçalo Pena et
Riccardo Puragliesi

13 Alternatives gratuites à Win-

dows: Ubuntu et Google
Laurent Kling
Résumé Helmholtz 3D, accompagnés de
résultats obtenus sur la nouvelle ar-
Cet article présente le paradigme de chitecture multicœur d’Intel: le Xeon 21 Nouveaux cours Java et .Net
programmation OpenMP pour les ma- 51XX Woodcrest. Notons enfin que Jean-Philippe Forestier
chines à mémoire partagée et propose les travaux présentés proviennent en
une méthodologie d’OpenMP-isation grande partie des travaux pratiques
des applications scientifiques sur les du cours postgrade High Performance 23 Programme des cours
machines SMP (Symetric Multi-Pro- Computing Methods de Ralf Gruber.
cessor) . Dans la première section, les
auteurs présentent en quelques mots le 28 Spécial été: IMAGES
contexte d’utilisation d’OpenMP ainsi Introduction Concours
qu’un petit historique. La deuxième de la meilleure nouvelle
section présente OpenMP dans les Il fut un temps, pas si lointain, où
grandes lignes. Une méthodologie l’utilisation de machines à mémoire
d’utilisation des machines à mé- partagée SMP n’était pas ouverte à tout
moire partagée avec OpenMP est un chacun: leur complexité matérielle Prochaines parutions
présentée dans la troisième section permettant une expression du paral- no délai parution
illustrée par un benchmarking sur le lélisme plus simple que les machines rédaction
simple exemple de la multiplication à mémoire distribuée avait un coût
matrice – matrice. Pour finir dans élevé. À cela il convient d’ajouter que 4 02.04.07 24.04.07
les cinquième et sixième sections, les toute l’énergie et le financement de la 5 03.05.07 22.05.07
6 07.06.07 26.06.07
auteurs présentent l’optimisation et recherche électronique ont été investis SP 28.06.07 28.08.07
la parallélisation avec OpenMP de dans la bataille de fréquence sur des 7 06.09.07 25.09.07
deux applications réelles: l’application processeurs mono-cœur que se sont 8 04.10.07 23.10.07
de physique des plasmas TERPSI- livrée les principaux fondeurs de micro- 9 01.11.07 20.11.07
CHORE, ainsi qu’un solveur de processeurs. 10 29.11.07 18.12.07

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

En novembre 1995, Intel annonçait le Pentium Pro, une Comment utiliser OpenMP ?
révolution dans l’évolution de son processeur phare capable
d’être monté sur une carte mère bi ou quadri processeurs. D’autres avant nous ont tenté d’établir une méthodolo-
Mais en partie à cause de la découverte du désormais célèbre gie d’OpenMP-isation (voir par exemple l’excellent cours de
flag erratum (un bug dans la conversion entier-float causant l’IDRIS [13]). Il nous semblait toutefois que la présentation
une erreur overflow), le chip n’a pas connu le succès qu’il aurait des résultats en terme de performances (orienté Calcul à
pu mériter, à l’exception de celui d’avoir été le précurseur Hautes Performances) en gardant le même exemple du début
de toute la série qui a suivi [1]. Aujourd’hui la guerre des à la fin, passant en revue les améliorations, les pièges à éviter,
GHz est terminée et on assiste à un revirement de stratégie etc. avait encore un sens.
des principaux fondeurs avec l’apparition des architectures Il existe deux cas d’utilisation: sur un code existant, ou
multi-processeurs à multi-cœurs, doublées de compilateurs from scratch. Dans le premier, on part d’un code séquentiel
qui s’améliorent de version en version, autorisant la pro- mono-thread qu’on parallélise au niveau des boucles en sui-
grammation à mémoire partagée à un coût relativement vant une méthodologie que nous allons voir plus loin; c’est
faible. OpenMP est l’un des paradigmes de programmation l’approche à grain fin (fine grain). Dans le second, on part dès
parallèle adapté aux machines de type SMP, c’est l’objet de le début en ayant en tête que la machine cible sera à mémoire
cet article. partagée; c’est l’approche à grain fort (coarse grain). Il est utile
Ces briques multi-processeurs multi-cœurs sont les nœuds de rappeler ici que les drapeaux de compilation (compilation
des clusters de type Beowulf actuels (clusters construits sur flags) disponibles sont nombreux et variés [9], nous en reste-
la base d’éléments grand public). Le 6ème ordinateur le plus rons à l’optimisation –O0 (sans optimisation), –O1, –O2 (les
puissant de la planète – et le premier Beowulf derrière les deux sont similaires sur les architectures IA-32 et EM64T et
spécifiques BlueGene d’IBM et le RedStorm de Cray Inc.– est correspondent à une optimisation avec inlining, propagation
un super-ordinateur de ce type [2]. Ainsi la programmation de constantes, allocation globale des registres, propagation
à plusieurs niveaux de parallélisme est maintenant à portée des attributs des routines, analyse des adresses des variables,
de bourse de tout un chacun; on utilise le paradigme de élimination des sections mortes, suppression des variables
programmation parallèle MPI entre les nœuds et OpenMP non référencées, optimisation inter routines au niveau d’un
à l’intérieur des nœuds. fichier, récursion de queue) et –O3 (optimisation maximale:
–O2 avec en plus une analyse plus agressive des dépendances,
modification des boucles, pré-extraction des données) ainsi
OpenMP, késako ? qu’à la vectorisation des boucles via les flags –xX où X dépend
de l’architecture cible (–xT sur un Woodcrest, –xW sur un
Le standard OpenMP (Open Multi-Processing [3]) est une Xeon, –xB sur un PentiumM, etc.).
interface de programmation (API) et a été publié pour la Afin de pouvoir mesurer les peformances, les effets des
première fois en 1997 pour le langage Fortran. Aujourd’hui, optimisations proposées, nous utiliserons le temps et le rap-
OpenMP est aussi disponible pour les langages de pro- port MFlops/s comme métriques. Ces deux valeurs peuvent
grammation C/C++. OpenMP consiste en un ensemble de se mesurer grâce à des outils spécifiques non intrusifs (Intel
directives de compilation, d’une bibliothèque de fonctions VTune [10] ou la simple routine gettimeofday() précise à
ainsi qu’un ensemble de variables d’environnement. Les la micro-seconde près) ou intrusifs (Parallel API [5]). Dans
performances obtenues avec OpenMP dépendent donc de cet exemple, on mesurera le temps à la microseconde près
deux facteurs: d’une part d’un choix de partage de données et on calculera le rapport MFlops/s avec la connaissance du
judicieux de la part du programmeur, d’autre part de la nombre d’opérations. Dans la suite de cet article nous allons
qualité du compilateur (interprétation des directives de prendre l’exemple simple de la multiplication matricielle (les
compilation). Tous les exemples décrits dans cet article ont éléments des matrices sont des nombres à virgule flottante
été compilés avec la dernière version du compilateur d’Intel double précision, 64 bits):
(9.1e) sur un nœud multi-cœurs du cluster Pleiades2 [4] C=A*B
(Intel Xeon 5150 Woodcrest, 2.66 GHz, 8 GB RAM), sauf La taille des matrices choisie est de 1200 x 1200 (division
lorsque l’auteur le précise explicitement. entière par 16 possible pour l’exemple). L’exemple sera implé-
Une machine à mémoire partagée se programme toujours menté en Fortran 771 et compilé à l’aide du compilateur Intel
de la même façon: plusieurs threads résidant sur des process Fortran 9.1e sur l’un des nœuds de Pleiades2+ (la nouvelle
elements différents possèdent un espace propre et accèdent extension multi-processeurs multicœurs du cluster Pleiades).
à un espace de données partagé; entraînant les problèmes L’exemple est implémenté suivant le listing 1: le listing 1
traditionnels d’aléas de données. L’utilisation d’OpenMP présente la séquence k – j – i, il suffit d’inverser les indices
ne change en rien cet état de fait et demande toujours du pour obtenir les résultats présentés dans cet article.
programmeur une connaissance parfaite du problème qu’il Tous les résultats sont donnés en forçant 4 threads sur le
veut résoudre en terme algorithmique ainsi que des bases de nœud, soit un thread par cœur
l’architecture cible. OpenMP peut améliorer sensiblement (export OMP–NUM–THREADS = 4).
la performance d’une application, mais peut malheureuse-
ment aussi, mal utilisé, la fusiller. C’est ce que nous allons
notamment voir dans la suite de cet article.

1
Ce choix se justifie par une visibilité plus claire de l’impact du choix de l’ordre des indices. Une analogie avec le langage C/C++ est aussi
plus aisée. La performance obtenue avec la syntaxe de base Fortran 90/95 sera présentée en fin de section

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

Première étape: nent la meilleure performance; la meilleure boucle est celle

compilation sans optimisation qui permet le meilleur agencement des données en mémoire.
La compilation sans optimisation (flag –O0) permet au Plus précisément, le principe du déroulement de boucle exige
programmeur de pointer les problèmes et surtout d’obser- de laisser en cache des grandeurs qui peuvent être utilisées
ver les améliorations qu’il peut faire lui-même à la main, plus d’une fois. Ainsi, on augmente le Va [6] (rapport entre
notamment l’ordre des indices dans les boucles. Le listing 2 le nombre d’opérations et le nombre de LOAD’s en mémoire
présente les résultats pour cette première étape. haute). Dans l’exemple qui suit, on déroule sur 8 niveaux,
cela signifie (voir listing 7 sans les directives OpenMP) que
n = 1200
c Initialisation de A, B et C les grandeurs b et c restent en cache (du fait de leur grandeur),
a est utilisé 8 fois, Va = 16 au lieu de Va = 2.
call initialise(n*n,A,1.0d0) La compilation se fait à l’aide du flag de compilation –O1.
call initialise(n*n,B,1.0d0) Le déroulement de boucle permet d’effectuer plus d’opéra-
call initialise(n*n,C,1.0d0)
tions avec les mêmes données en cache. Le programmeur
c Depart du chrono doit par contre toujours être clair quant aux indices qu’il
ti = second() parcourt.
c Multiplication C=A*B
c en suivant la sequence kji t Mflop/s
do k=1,n
kji 0.3889E+01 888.6
do j=1,n
do i=1,n kji2 0.2133E+01 1620.0
C(j,k)=C(j,k)+A(j,i)*B(i,k) kji4 0.1596E+01 2166.0
enddo
enddo kji8 0.1584E+01 2182.0
enddo kji16 0.2077E+01 1664.0
c Arret du chrono
ti = second ()-ti listing 3 – Résultats du programme MMIJK avec le flag de
mflops_kji = 2.0*(n/100.)**3/ti
compilation -O1 et un unroll manuel
listing 1 – Multiplication matrice*matrice
Le listing 3 démontre qu’un unroll manuel améliore les
t Mflop/s
performances de la meilleure des boucles. Ainsi on gagne en-
core 1.3 GFlops/s sur le meilleur arrangement des indices.
kji 0.3849E+02 89.78
kij 0.2339E+02 147.7
ikj 0.2367E+02 146.0
Troisième étape:
vectorisation des boucles
jki 0.2927E+02 118.1
La documentation d’Intel [9] parle de innermost loops
ijk 0.4767E+02 72.50 parallelization c’est-à-dire de parallélisation sur la boucle
jik 0.4761E+02 72.59 intérieure alors que la parallélisation multi-threads s’appelle
DGEMM 0.1319E+00 29210.0 outer-most loops parallelization (voir sections 3.5.1 et 3.5.2).
L’option de compilation est –xX où X correspond à l’archi-
listing 2 – Résultats du programme MMIJK avec le flag de tecture cible. Cette optimisation permet d’utiliser au mieux
compilation–O0 les pipelines des processeurs en faisant plus d’opérations avec
les mêmes opérandes à chaque itération des boucles. On
Le listing 2 démontre que l’ordre des indices dans une n’utilisera cette option qu’après avoir effectué l’unrolling à la
boucle a une importance capitale dans la performance d’un main (étape 2). La performance maximale (2.18 GFlops/s) est
code. Les performances peuvent doubler. Ici le pire cas est de toujours celle obtenue grâce au unroll à 8 niveaux en regard
89 MFlops/s, le meilleur à 147 MFlops/s. Cela est dû au fait de la moins bonne (ordre j – i – k), de 354 MFlops/s, c’est
que les données restent en cache au lieu d’être à chaque fois 1.82 GFlops/s de plus donc un facteur de 5.
ramenées de la mémoire haute aux registres. DGEMM, la Notez que le compilateur est capable de faire un unroll sur
routine BLAS, élément de LAPACK[12] et écrite en assem- les boucles j–k–i et j–i–k mais n’arrive jamais à la performance
bleur, provenant de la bibliothèque MKL sur les machines obtenue grâce à un déroulement manuel. Les performances
Intel (ACML sur les machines à base d’AMD Opteron) a été obtenues en commençant par l’indice i sont de l’ordre de
mise comme référence. La performance de 29.21 GFlops/s 1 GFlops/s (soit moins de 40% de la performance avec un
est obtenue grâce aux instructions SSE2 et SSE3 (utilisant déroulement manuel), celles obtenues en commençant par
un chip spécifique permettant 4 résultats par temps de cycle l’indice k tournent autour de 600 MFlops/s (soit moins de
contre 2 normalement). Les plus perspicaces auront noté que 25% de la performance avec un déroulement manuel).
29.21 GFlops/s c’est 8 GFlops/s de plus que la performance Finalement, qu’on compile avec le flag le plus (–O3) ou
maximale du nœud (R∞ = 21.28 GFlops/s)2. le moins (–O1) agressif, seules les boucles commençant par
l’indice j gagnent en performance. Ceci se vérifie simplement
Seconde étape: entre les listings 5 et 4. Cela démontre que le compilateur
unroll manuel est capable de reconnaître certaines constructions et de les
Attaquons le déroulement de boucles (unrolling) manuel optimiser, malheureusement pas toutes.
de la meilleure des boucles k – i – j, celle dont les indices don-
2
R∞ = freq * NbrOperations per cycle * Nbrprocess elements = 2.66 * 109 * 2 * 4 = 21.28 GFlops/s

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

t Mflop/s Parallélisation manuelle avec OpenMP

kji 0.3543E+01 975.6 Le listing 7 présente le code source pour la boucle
kij 0.3543E+01 975.4 k–j–i avec un unroll manuel à 8 niveaux. Les directives de
ikj 0.5978E+01 578.2 programmation OpenMP sont placées à l’extérieur de la
jki 0.9751E+01 354.4 boucle à paralléliser. La boucle qui sera parallélisée est celle
ijk 0.5966E+01 579.3
qui suit immédiatement la directive !$OMP DO. La directive
SCHEDULE(GUIDED, 1) permet de choisir la façon dont
jik 0.2056E+02 168.1
seront agencées chacune des sections traitées en parallèle.
listing 4 – Résultats du programme MMIJK avec les flags de Dans ce cas précis, nous avons choisi un ordonnancement
compilation -O1 -xT sans unroll GUIDED i.e. que chaque thread reçoit une quantité de travail
balancé (load balancing) selon l’état de ce qui reste à faire.
t Mflop/s
L’auteur propose au lecteur intéressé de consulter les spécifi-
cations du standard OpenMP [11] pour plus de détails.
kji 0.3549E+01 973.7
kij 0.3548E+01 974.2
n = 1200
ikj 0.5988E+01 577.1 c Initialisation de A, B e t C
jki 0.1656E+01 2086.0
call initialise (n*n,A,1.0d0)
ijk 0.5987E+01 577.3
call initialise (n*n,B,1.0d0)
jik 0.1659E+01 2084.0 call initialise (n*n,C,1.0d0)

listing 5 – Résultats du programme MMIJK avec les flags de c Depart du chrono

compilation -O3 -xT sans unroll manuel ti = second ()

c Multiplication C = A * B
Quatrième étape: c en suivant la sequence kji et unroll de 8 niveaux

vectorisation et optimisation !$OMP PARALLEL DEFAULT(SHARED) PRIVATE( i , j )

!$OMP DO SCHEDULE(GUIDED, 1)
t Mflop/s do k=1,n , 8
do i =1,n
kji 0.1296E+01 2667.0
do j =1,n
kji2 0.1628E+01 2123.0 c(j,k)=c(j,k)+a(j,i)–b(i,k)
kji4 0.1493E+01 2314.0 c(j,k+1)=c(j,k+1)+a(j,i)–b(i,k+1)
c(j,k+2)=c(j,k+2)+a(j,i)–b(i,k+2)
kji8 0.1567E+01 2205.0 c(j,k+3)=c(j,k+3)+a(j,i)–b(i,k+3)
kji16 0.2077E+01 1664.0 c(j,k+4)=c(j,k+4)+a(j,i)–b(i,k+4)
c(j,k+5)=c(j,k+5)+a(j,i)–b(i,k+5)
listing 6 – Résultats du programme MMIJK avec les flags de c(j,k+6)=c(j,k+6)+a(j,i)–b(i,k+6)
c(j,k+7)=c(j,k+7)+a(j,i)–b(i,k+7)
compilation -O3 -xT et unroll manuel enddo
enddo
Le listing 6 présente le résultat après la compilation la enddo
!$OMP END DO
plus agressive et la vectorisation. La meilleure performance !$OMP END PARALLEL
obtenue est de 2.66 GFlops/s pour la boucle k–i–j ; ce qui
représente 50 % de R∞ ! Il est intéressant de voir ce que c Ar ret du chrono
donnent les performances d’un code sans aide au compilateur. ti = second () – ti
mflops_kji8 = 2 . 0*(n/100.)**3/ti
C’est ce qui est présenté au listing 5.
listing 7 – Multiplication matrice*matrice parallélisée avec
Cinquième étape: OpenMP
parallélisation
Une fois que le code est optimisé sur un processeur Le listing 8 démontre que la parallélisation par OpenMP
(ou un cœur), on peut commencer à paralléliser le code au sur 4 threads suit un speedup parfaitement linéaire (voir la
niveau du nœud (puisque nous sommes sur une machine comparaison avec les listings 2 et 3). Pour certaines boucles
SMP). Deux possibilités s’offrent au programmeur: laisser (i–k–j par exemple), on observe un speedup superlinéaire3.
le compilateur découvrir les sections qui peuvent être exé-
cutées en parallèle de façon sûre (safely parallelized sections), Auto-parallélisation
c’est l’auto-parallélisation, ou alors choisir manuellement la Cette étape consiste à laisser le compilateur chercher (et
meilleure façon de partager le travail entre chaque thread, trouver) les boucles qui peuvent s’exécuter de façon indé-
c’est la parallélisation avec OpenMP. Nous allons démon- pendante (embarrassingly parallel). Pour cela, on invoque le
trer que cette seconde façon est la meilleure, même pour le compilateur à l’aide du flag de compilation –parallel. Cette
programme le plus simple, puisqu’on n’est jamais mieux servi opération se fait après l’optimisation –O2 (ou –O3) A noter que
que par soi-même. si le programmeur a lui-même procédé à la parallélisation de
ses boucles en utilisant OpenMP (comme décrit à la section

3
Cette dernière remarque permet aussi à l’auteur de mettre en garde le lecteur ne jurant que par le speedup: si la version monoprocesseur
affiche une performance médiocre, il est relativement aisé d’obtenir ce genre d’effet !

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

précédente), l’invocation de -parallel -openmp ne change Sixième étape:

pas le programme final: le compilateur parallélise d’abord les optimisation du code parallèle
boucles OpenMP, puis autoparallélise le reste. L’autoparallé- Que se passe-t-il maintenant si on force le compilateur à
lisation ne peut se faire qu’avec le flag d’optimisation –O1 au optimiser les boucles parallélisées avec OpenMP ?
minimum, la comparaison entre la parallélisation OpenMP et Le listing 12 présente les résultats du code parallélisé
le niveau –O0 est donc difficile. On remarque tout de même avec OpenMP et le flag d’optimisation –O3. La performance
que le maximum que puisse faire le compilateur sans unroll maximale obtenue est de 7.73 GFlops/s pour la séquence
manuel est 2.84 GFlops/s (voir listing 9) alors que si le pro- k–i–j avec un unrolling de 8. Les résultats sont exactement
grammeur a suivi les étapes proposées par l’auteur aux points les mêmes qu’avec le flag de compilation –O2 !
précédents, notamment l’unrolling manuel, le compilateur Cela signifie que le compilateur n’est plus capable de
arrive à des performances de l’ordre de 7.77 GFlops/s pour un faire mieux. Finalement, le listing 13 présente les résultats
unrolling manuel à 8 niveaux. L’autoparallélisation optimisée du code parallélisé avec OpenMP, le flag d’optimisation
de la façon la plus agressive donne des résultats proches de ceux –O3 et la vectorisation des boucles intérieures avec –xT. La
d’OpenMP optimisé pour certains indices de boucle. Voir le performance maximale obtenue est de 8.04 GFlops/s pour
listing 11. Il est intéressant de voir que la suite d’index la plus la séquence k–i–j avec un unrolling de 8. C’est le meilleur
utilisée i–j–k donne la performance la plus médiocre. résultat possible4.
t Mflop/s
t Mflop/s
kji 0.9550E+01 361.9
kji 0.2847E+01 1214.0
kij 0.5862E+01 589.6
kij 0.1211E+01 2853.0
ikj 0.6714E+01 514.7
ikj 0.3092E+01 1118.0
jki 0.7697E+01 449.0
jki 0.2470E+01 1399.0
ijk 0.1531E+02 225.7
ijk 0.8818E+01 391.9
jik 0.1396E+02 247.5
jik 0.7624E+01 453.3
kji2 0.5589E+01 618.3
kji2 0.6555E+00 5273.0
kji4 0.5504E+01 628.0
kji4 0.5216E+00 6626.0
kji8 0.5800E+01 595.9
kji8 0.4468E+00 7736.0
kji16 0.5602E+01 616.9
kji16 0.5124E+00 6745.0
listing 8 – Résultats du programme MMIJK avec le flag de
compilation -O0 -openmp listing 12 – Résultats du programme MMIJK avec le flag de
compilation -O3 -openmp
t Mflop/s
kji 0.2871E+01 1204.0 t Mflop/s
kij 0.1213E+01 2848.0 kji 0.1166E+01 2964.0
ikj 0.2409E+01 1435.0 kij 0.1165E+01 2967.0
jki 0.2471E+01 1398.0 ikj 0.3190E+01 1083.0
ijk 0.6027E+01 573.5 jki 0.6315E+00 5473.0
jik 0.5439E+01 635.4 ijk 0.6405E+01 539.6

listing 9 – Résultats du programme MMIJK avec le flag de jik 0.6076E+00 5688.0

compilation -O1 -parallel sans unroll manuel kji2 0.6163E+00 5607.0
t Mflop/s kji4 0.5256E+00 6576.0
kji2 0.6550E+00 5276.0 kji8 0.4298E+00 8042.0
kji4 0.5193E+00 6655.0 kji16 0.4658E+00 7420.0
kji8 0.4448E+00 7769.0 listing 13 – Résultats du programme MMIJK avec le flag de
kji16 0.5081E+00 6802.0 compilation -O3 -openmp -xT
listing 10 – Résultats du programme MMIJK avec le flag de
compilation -O1 -parallel avec unroll manuel Septième étape:
toutes voiles dehors
t Mflop/s Pour finir, regardons ci-après ce que donne le résultat
kji 0.1177E+01 2936.0 compilé avec –O3–xT–parallel– openmp. Le compilateur
kij 0.1171E+01 2952.0 unroll parallélise les boucles selon OpenMP, puis auto-pa-
ikj 0.2449E+01 1411.0 rallélise ce qui reste et finalement vectorise.
jki 0.4249E+00 8134.0 Le listing 14 présente les mêmes performances que celles
ijk 0.3780E+01 914.4 obtenues avec le code parallélisé à la main avec OpenMP
jik 0.4249E+00 8133.0
(listing 13). L’auto-parallélisation sur un code optimisé à la
main n’améliore pas les performances même si le compilateur
listing 11 – Résultats du programme MMIJK avec le flag de annonce (lors de la compilation) qu’il a réussi à auto-paral-
compilation -O3 -xT -parallel sans unroll manuel léliser des boucles partielles.
4
Il s’agit aussi tout simplement de la performance maximale R1 du processeur vectoriel NEC SX-5 (R1 = 8 GFlops/s)
FI 3 – 27 mars 2007 – page 
Optimisation et parallélisme avec OpenMP: l'approche grain fin

t Mflop/s Moralité
kji 0.1166E+01 2964.0 Les résultats présentés ci-dessus démontrent – à l’aide
kij 0.1166E+01 2964.0 d’un code source extrêmement simple – qu’en suivant une
ikj 0.2981E+01 1159.0 stratégie d’optimisation on peut arriver à des performances
jki 0.6556E+00 5272.0 de l’ordre de 37 % de la performance maximale R∞ du pro-
ijk 0.6453E+01 535.5 cesseur considéré. Sans optimisation manuelle, en plaçant
jik 0.5951E+00 5807.0
une totale confiance dans le compilateur, on obtient un
résultat proche (flags de compilation -parallel -xT -O3)
kji2 0.6160E+00 5611.0
8.1 GFlops/s au mieux dans le cas où l’utilisateur a aidé
kji4 0.5263E+00 6566.0
manuellement le compilateur (1.4 GFlops/s au pire !). Le
kji8 0.4300E+00 8037.0
compilateur est capable de reconnaître certaines boucles et de
kji16 0.4658E+00 7419.0
changer l’ordre des indices. Malheureusement pas partout.
listing 14 – Résultats du programme MMIJK avec le flag de Les performances de l’auto-parallélisation démontrent
compilation -O3 -xT -parallel-openmp que le compilateur est capable de reconnaître certaines bou-
cles et de les paralléliser exactement comme on le fait avec
Et Fortran 90 alors ? OpenMP. Pour certains ordres d’indice, le compilateur n’en
Afin d’être complet, donnons quelques performances est pas capable. Il est aussi à noter que suivant la complexité
pour qui utiliserait les nouvelles fonctionnalités offertes par du calcul à l’intérieur de la boucle externe, le compilateur ne
Fortran 90. La principale nouveauté apparue avec Fortran sera pas capable de l’auto-paralléliser. C’est ce que nous allons
90 sur Fortran 77 était la gestion simplifiée des tableaux et voir dans la dernière partie de cet article avec l’application
l’apparition de nouvelles fonctions intrinsèques. Dans le cas de physique des plasmas TERPSICHORE d’abord, avec un
de notre petit exemple, la multiplication matricielle se note solveur de Helmholtz 3D ensuite.
tel que présenté au listing 15. Cela montre surtout que si le programmeur a une
connaissance exacte de l’algorithmique de son code, il peut
n = 1200 alors parfaitement maîtriser la localisation de ses données.
! Initialisation de A, B e t C
call initialise(n*n,A,1.0d0)
call initialise(n*n,B,1.0d0)
call initialise(n*n,C,1.0d0) Premier exemple: TERPSICHORE
ti = second() TERPSICHORE est un code MHD développé au CRPP
C = matmul(A,B)
ti = second ( ) – t i dès 1988 [7].
mflops–rate = 2.0*(n/100.)**3/ti Kinetic theory

Plasma
listing 15 – Multiplication matricielle en Fortran 90 physics Two-fluid model
resistive MHD Equilibrium

MHD ideal MHD Stability

Il est évident que cette notation à la Matlab est bien plus
lisible pour un humain, mais qu’en est-il des performances ? Fig. 1 – TERPSICHORE étudie la stabilité MHD idéale des
Elles sont catastrophiques: appareils de fusion thermonucléaire
t Mflop/s
F90 0.1096E+02 315.3 La physique de TERPSICHORE
DGEMM 0.1203E+00 29120.0 TERPSICHORE est un code MHD idéal (MagnetoHy-
droDynamics), nommé ainsi en référence à la déesse grecque
listing 16 – Performance MATMULT en Fortran 90 avec flags
de la danse et de la poésie. Il recherche la stabilité des géo-
-O3 -xT
métries tridimensionnelles à l’aide d’un équilibre provenant
d’un autre code nommé VMEC [8]. L’équilibre est trouvé
On a laissé la performance de DGEMM de côté. Ce n’est
en résolvant le système non linéaire:
pas moins d’un facteur 100 qui différencie les deux métho-
des! Une OpenMP-isation devrait améliorer la performance p = j × B,       j = × B,        · B = 0,
(voir listing 17) avec p, j, B: la pression, la densité de courant et le champ
magnétique respectivement. Ensuite, pour déterminer la
t Mflop/s stabilité du plasma, les équations linéarisées MHD sont
F90// 0.2905E+01 1190.0 traitées par TERPSICHORE. Elles peuvent être écrites sous
la forme variationnelle :
listing 17 – Performance MATMULT en Fortran 90 avec flags
-O3 -xT - openmp dWp + dWv − w2 dWk = 0,
avec dWp l’énergie potentielle du plasma, dWv l’énergie ma-
Et Fortran 95 alors ? gnétique du vide autour du plasma, dWk l’énergie cinétique
Le standard Fortran 95 a introduit des nouvelles fonction- et w2 la valeur propre du système. Si w2 < 0 alors la configu-
nalités parallèles. C’est le cas notamment du bloc FORALL. ration est instable. La fig. 2 montre un exemple d’équilibre
Dans le cas de la multiplication matricielle, cela n’a pas le traité avec TERPSICHORE.
moindre sens.

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

Organisation du code
La structure de TERPSICHORE est organisée en six
principales routines: eqinvm, veqrec, mtaskb, sta-
bin, mtaska et mtasks–ap. Les deux premières établissent
l’interface avec la sortie du code VMEC et reconstruisent
les variables d’équilibre MHD. La routine mtaskb est
responsable principalement de la mise en correspondance
du système de coordonnées de VMEC et du système de
coordonnées de Boozer ainsi que du calcul de la métrique.
La construction des matrices de stabilité est effectuée dans
stabin. Le solveur PAMERA calcule les valeurs propres w2
et les modes. Ce solveur prend peu de temps (ne domine
pas) pour le petit cas test considéré.

Fig. 2 – un exemple de configuration: le W7-X à Greifswald,

Allemagne. Performance de TERPSICHORE
Le code a été optimisé dans le passé pour les machines
Une des spécialités de TERPSICHORE est son système vectorielles, spécialement pour le NEC SX-5 dans sa dernière
de coordonnées. La MHD idéale prédit que les lignes de version et donc, la structure montrée à la fig. 3 est plus pré-
champ magnétique sont sur des surfaces à pression constante, cisément donnée par:
ce qu’on appelle les surfaces de flux, et que ces lignes de
champ magnétique sont gelées dans le plasma, c’est-à-dire do s
do theta
que si les particules bougent, les lignes de champ magnétique do psi
bougent avec elles. Avec cette connaissance, il est possible variable(s,theta,psi)
d’introduire des coordonnées du flux magnétique, ou coor- end do
end do
données de Boozer. Ces coordonnées de Boozer utilisent end do
une étiquette du flux s [0, 1] comme variable radiale et les
variables angulaires polaires et toroïdales q, y de façon à ce La performance sur la machine NEC SX-5 pour le run
que les lignes de champ deviennent droites dans le système de test choisi (Large Helical Device) est de 30s. Cette perfor-
de coordonnées choisi. Si l’on veut ensuite considérer une mance devient l’objectif à atteindre en optimisant le code sur
particule dans une géométrie donnée, on sait d’avance qu’elle l’architecture Xeon 5100 (Woodcrest) de Pleiades2.
ne va pas bouger le long de la variable radiale s du fait de sa Le tab. 1 présente le profile de l’application non opti-
condition de gel. Il s’agit là de l’un des nombreux avantages misée sur une architecture XEON. Quatre routines attirent
de ce nouveau système de coordonnées. immédiatement l’attention : lgikvm, lamcal, metric et
fourin. La suite de cette section concernera l’optimisation
de ces quatre routines principalement.

main routine subroutine subroutine time

eqinvm 0.12
do s
do theta veqrec cospol
do psi
0.64
... lgikvm 35.45
end do
mtaskl lamcal 24.78
... s
end do mtaskb vmtobo
... 24.16
end do extint 5.73
metric 15.73
bophys 4.96
stabin 0.14
mtaska fourin 10.01
Fig. 3 – structure schématique
Total 132.97
Le code
L’utilisation de coordonnées magnétiques permet de cal- Tab. 1 – TERPSICHORE: profile obtenu sur un Xeon 64 HT
culer chaque surface de flux (c’est-à-dire une position radiale) (Pleiades2) de l’application non optimisée
de façon indépendante. Ainsi, la structure générale est très
similaire à un code dans lequel il y aurait des boucles sur Optimisation
chaque variable radiale et ce, sur chaque sous-routine. Sché- Performance sur le Xeon 5160 (WoodCrest)
matiquement, cette structure est présentée dans la fig. 3. Tout d’abord, nous avons voulu savoir où commencer
notre optimisation. La fig. 4 ainsi que le tab. 2 présentent la

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

performance du code obtenue sur un Xeon 5160 (Woodcrest). Notons tout d’abord que toutes les exécutions montrent
On voit facilement que le compilateur utilisé (ifort version une meilleure performance que celle obtenue sans optimisa-
9.1e) peut améliorer les résultats d’environ 50%, mais, tion (valeur initiale), c’est-à-dire que l’optimisation à la main
comme nous l’avons vu plus haut, il est recommandé d’aider avec le flag de compilation (–O1 par exemple) était meilleure
le compilateur manuellement plutôt que d’accorder une trop d’environ 15% que ce dont est capable le compilateur (avec
grande confiance en lui. Dans tous les cas, la nature du code –O3 –xT). Deuxièmement, on observe qu’il n’y a quasiment
laisse présager d’excellentes optimisations possibles. pas de différence de performance sur la version optimisée
manuellement en changeant de flag de compilation (excepté
avec –parallel). Ceci nous encourage donc, comme précisé
dans la première partie de cet article, à paralléliser le code ma-
nuellement pour voir si on peut encore aider le compilateur
et obtenir les meilleures performances possibles (meilleures
que ce que pourrait faire le compilateur seul).

Fig. 4 – performance sur un nœud Xeon 5160 (WoodCrest)

avant l’optimisation

-O1 133s
-O3 124s
-O3-xT 96s Fig. 6 – la performance après optimisation sur un processeur
-O3 -xT -parallel 65s au niveau –O1
Tab. 2 – TERPSICHORE: performance sur un nœud Xeon
opt
5160 (WoodCrest) avant l’optimisation
-O1 133s 85s
-O3 124s 84s
Travailler sur l’optimisation -O1 -O3 -xT 96s 84s
Comme déjà vu plus haut, il y a quelques changements
-O3 -xT -parallel 65s 51s
qui doivent être apportés manuellement à ce niveau d’opti-
misation. Alors que le code a été déjà hautement optimisé, Tab. 3 – TERPSICHORE: performance après optimisation
mais pour une machine vectorielle, le gain principal que
l’on peut attendre vient d’une simple inversion des indices Stratégie de parallélisation
des boucles DO. La fig. 5 présente un exemple d’une telle Le code sera parallélisé en utilisant les directives OpenMP
inversion. D’autres possibilités ont été utilisées comme la sur le code optimisé sur un processeur.
division des grandes boucles, le déroulement dans des plus Nous avons observé que, pour ce code, la meilleure
petites boucles, etc. mais l’inversion des indices reste le plus stratégie de partage des données entre les threads OpenMP
efficace. Les résultats après cette première optimisation sont était (dans l’ordre):
présentés à la fig. 6 et au tab. 3. z !$OMP DO sur la boucle la plus exterieure;

Fig. 5 – La principale accélération (speedup) a été obtenue par inversion des indices dans les boucles DO

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

Fig. 7 – un exemple d’utilisation de la directive de compilation !$OMP SECTIONS dans une boucle indépendante (embarrassingly
parallel loop)
z division des grandes boucles indépendantes et placement Résultats
dans des directives !$OMP SECTIONS; Comme petit résumé de ce qu’on a pu observer de la
z combinaisons de petites boucles; parallélisation avec OpenMP, nous pouvons établir que –>
z appels OpenMP imbriqués; l’optimisation et la parallélisation à la main est meilleure
La fig. 7 montre un exemple de partage en sections que le compilateur –> sur des codes instables peut être aussi
indépendantes d’une boucle et la performance associée (en rapide que sur des codes stables avec OpenMP –> le Xeon
secondes). 5160 (WoodCrest) est comparable au NEC SX-5... mais il
Après la parallélisation avec OpenMP, la performance a est 30 fois moins cher. Le tab. 5 présente un petit résumé des
augmenté d’un facteur de 3 en comparaison de la version performances obtenues selon les flags de compilation ainsi
optimisée non parallèle (fig. 8 et tab. 4). que le paradigme de parallélisation utilisé (autoparallélisation
On peut observer qu’une parallélisation manuelle est ou manuellement avec les directives OpenMP).
meilleure que ce que peut faire le compilateur : -O1 -openmp
est environ 2 fois plus rapide que -O3 -xT -parallel sur
main routine routine routine non opt sp5 opt
le code non optimisé. Notons aussi qu’il n’y a quasiment
eqinvm 0.12 0.13
pas de différence entre les flags de compilation et que la
performance ainsi obtenue est comparable à celle obtenue veqrec cospol 0.64 0.0.64
sur le NEC SX-5. lgikvm 35.45 13.15
mtaskl lamcal 24.78 17.77
mtaskb vmtobo 24.16 19.35
extint 5.73 5.78
metric 15.73 6.18
bophys 4.96 5.50
stabin 0.14 0.14
mtaska fourin 10.01 7.74
Total 132.97 84.69
Tab. 5 – résumé des performances des quatre principales
routines posant problème

Il est intéressant de noter que le gain de performance

Fig. 8 – performance après la parallélisation avec OpenMP entre la version OpenMP-isée et la version optimisée est d’un
opt OMP facteur de 3.1 pour 4 threads OpenMP, signifiant ainsi que
-O1 133s 85s 31s la scalabilité est bonne. Difficile de faire mieux.
-O3 124s 84s 31s
-O3 -xT 96s 84s 28s
-O3 -xT -parallel 65s 51s 27s Second exemple: Solveur de
NEC SX-5 30s Helmholtz 3D
Tab. 4 – performance après parallélisation. Les pourcentages La problématique
sont relatifs à avant la parallélisation avec le même flag de L’équation de Helmholtz peut être utilisée pour résoudre
compilation des problèmes complexes par exemple un flux thermal induit
par les forces de flottaison (buoyancy forces) à l’intérieur d’un

FI 3 – 27 mars 2007 – page 

Optimisation et parallélisme avec OpenMP: l'approche grain fin

volume fermé où la différence de température est imposée Hi Li Hi–1 = L1* (8)

entre les bords. et il est possible de démontrer que Hi et Li sont proportion-
En partant des équations de Navier-Stockes et de l’équa- nels à la matrice identité dans les directions j ≠ i.
tion de l’énergie sous l’hypothèse de Boussinesq, il est possible
de diviser le problème initial comme suit [14, 15]: Optimisation sur un processeur
[Première étape] Avant de commencer une quelconque parallélisation du
(1) code, il est très important, comme nous l’avons vu précédem-
ment, de tirer le maximum de performance de l’application
(2) sur un processeur en optimisant manuellement. L’objectif
étant d’obtenir un exécutable qui offre une performance
Du
a* • n = $u • n in D7 similaire, quels que soient les flags de compilation, même
Dt (3)
si le déroulement des boucles et la vectorisation des boucles
[Seconde étape] internes sont faits par le compilateur.
Du Un bon point de départ est d’identifier la routine qui
$u = a* in 7 prend le plus de temps sur l’entier du code, parce qu’une
Dt
(4) infime amélioration de cette routine entraîne une amélio-
plus les conditions de bord. ration globale non négligeable. En utilisant l’outil gprof 6
La méthode de projection-diffusion brièvement décrite couplé avec l’option de compilation -qp(-pg), on obtient
ci-dessus entraîne quatre problèmes de Helmholtz utiles pour directement un profile du code. Le tableau 6 présente le
résoudre chaque itération: un pour chaque composante carté- profiling de l’entier du solveur de Navier-Stokes.
sienne de la vitesse et un pour la température. Puisque chaque À première vue, il semble que la routine helm3col (le
variable primitive est étendue en polynôme de Chebyshev, solveur de Helmholtz) est la partie du code la plus gour-
pour le cas présent, une méthodologie pseudo-spectrale est mande en temps. C’est sur cette routine que nous allons
utilisée en imposant une erreur nulle aux points de la grille nous concentrer.
(plus précisément aux points de Gauss-Lobatto), une grille Le tab. 7 présente les performances en terme de temps
de 1693 points est suffisante pour exécuter une simulation passé pour une seule itération enregistrée sur un processeur
numérique directe (DNS Direct Numerical Simulation) dans de la machine Pleiades2 (Xeon 64, 2.8 GHz). La première
un régime turbulent. colonne présente le cas d’une compilation non agressive -O1
Donc, la matrice discrète du problème pour une quantité sans optimisation à la main, la seconde présente le résultat
scalaire générale peut être écrite comme: obtenu avec une compilation agressive -O3 -xW, ensuite
*2
Øn+1 – a(Dt)Øn+1 = f * in W (5) les résultats présentent les mêmes flags de compilation mais
ou, en séparant l’opérateur Laplacien discret modifié avec une optimisation à la main. Par optimisation manuelle,
(L*x + L*y + L*z – a(Dt)I) Øn+1 = F* (6) nous entendons simplement que toutes les boucles DO ont été
En outre, résolvant un problème de valeur propre – vec- modifiées afin d’aider le compilateur à les traiter (réordon-
teur propre de l’opérateur Laplacien discret modifié pendant nancement des index, déroulement, division, etc.).
une étape de prétraitement, alors le problème peut être résolu
simplement en utilisant la multiplication matrice – matrice NH-NACO NH-ACO NH-ACO NH-ACO
pour chaque direction cartésienne. En fait, le problème [s/iter] [s/iter] [s/iter] [s/iter]
final établit: 33.39(33.32) 6.90(4.42) 3.60(3.63) 3.10(3.12)
HxHyHz [Lx + Ly + Lz + a(Dt)I] Hx–1Hy–1Hz–1Øn+1 = F1*
Tab. 7 – Différentes performances obtenues par modification
(7)
des flags de compilation et optimisation manuelle
où

Each sample counts as 0.01 seconds

% cumulative self calls self total name
time seconds seconds s/call s/call
38.20 447.73 447.73 200 2.24 2.82 helm3col–
6.53 524.26 76.53 mkl–blas–p4–dgcopyan
6.36 598.79 74.53 1252 0.06 0.06 permyz–
6.17 671.07 72.28 1252 0.06 0.06 permyz–
5.28 732.91 61.84 mkl–blas–p4–dcopy
5.00 791.51 58.60 mkl–blas–p4–dgcopybn
...
Tab. 6 – Profile plat en utilisant gprof avec les options de compilation adéquates. La routine qui prend le plus de temps est le solveur
lui-même: helm3col
6
sur une architecture sur laquelle on peut faire confiance à cet outil (ndt)

FI 3 – 27 mars 2007 – page 10

Optimisation et parallélisme avec OpenMP: l'approche grain fin

En outre, nous avons testé la capacité HT (HyperThrea- Ainsi il est fortement recommandé de réécrire les sous-rou-
ding) (voir [9]) des processeurs Xeon afin d’obtenir un gain tines BLAS 1 (dscal, dcopy, etc.) en utilisant les directives
pour l’exécution d’une instruction parallèle très simple (la ca- OpenMP. Le tab. 10 présente la loi d’échelle pour la routine
pacité HT autorise une augmentation infime de performance dcopy (de la bibliothèque MKL version 8.0) ainsi que pour
sur un Pentium 4). Le tab. 8 résume ce qui a été obtenu en celle écrite à la main ucopy sur 1693 éléments de matrice.
utilisant le flag de compilation –parallel (l’autoparalléli- Cela même si les BLAS sont implémentés dans MKL avec
sation automatique du compilateur). Les deux versions du les directives OpenMP.
compilateur sont à nouveau comparées.
compiling 1cpu 2cpus 4cpus
NH-ACO NH-ACO options [s/iter] [s/iter] [s/iter]
[s/iter] [s/iter] -parallel 2.07 1.37 1.07
33.00(4.85) 3.15(3.21) -openmp 2.07 1.20 0.87
Tab. 8 – Influence de l’HyperThreading sur les performances
avec l’auto-parallélisation du compilateur Tab. 9 – Performance auto-parallélisation vs. directives ma-
nuelles OpenMP sur 1, 2 et 4 cœurs d’un nœud multi-CPU’s
Il apparaît clairement que la toute dernière version du multi-cœurs Woodcrest avec le code totalement parallélisé
compilateur Intel (la version 9.1e) n’est pas capable de corri-
ger les erreurs d’écriture du code, les performances sont ainsi
subroutine 1cpu [s] 2cpus [s] 4cpus [s]
dégradées drastiquement. En effet, dans ce cas, l’auto-paral-
lélisation semble annuler toutes les autres optimisations et dcopy 2.56 x 10 –2
2.56 x 10 –2
2.56 x 10–2
donne une performance obtenue avec une optimisation basse ucopy 2.54 x 10–2 1.61 x 10–2 1.58 x 10–2
–O2. D’un autre côté, la version précédente du compilateur
(la version 9.0 disponible sur Pleiades2) donne des résultats Tab. 10 – Loi d’échelle des deux routines dcopy (de MKL) et
bien meilleurs même s’ils ne sont pas comparables avec une ucopy (écrite avec OpenMP)
optimisation manuelle. Notons finalement que la compila-
tion du code optimisé manuellement avec l’une ou l’autre des
versions du compilateur donne les mêmes résultats. Quelques observations
Parallélisation Dans cette dernière section, nous donnons quelques
L’étape suivante est la parallélisation du code sur une conseils sur la base d’observations faites lors de l’optimisation
machine SMP. Pour ce faire, nous utiliserons les nouvelles et la parallélisation d’applications à l’aide d’OpenMP.
architectures multi-cœurs d’Intel Woodcrest. Tout ce qui va De façon générale, la procédure proposée est celle com-
suivre a été obtenu sur un seul nœud du calculateur Pleia- munément utilisée en génie logiciel lorsqu’il s’agit d’optimiser
des2+. Les résultats qui vont suivre présentent uniquement une application: Monitoring –> Optimisation –> Tests.
une parallélisation manuelle; l’avantage de l’utilisation des Dans cette optique, il faut disposer d’une bonne mesure
directives OpenMP à la place de l’auto-parallélisation va du temps. Des outils existent (VTUNE, gprof, etc.). Mais
être démontrée. nous avons pu observer que le meilleur monitoring reste
L’utilisation des directives OpenMP a permis de paralléli- l’utilisation de la fonction gettimeofday(). Le profiling
ser des opérations complexes qui concernent principalement gprof donne des résultats erronnés.
une quadruple boucle imbriquée où une multiplication ma- La version du compilateur ainsi que celle des versions
tice – matrice avec des indices shiftés est effectuée entre une utilisées peuvent jouer un rôle. Dans cette optique, les ingé-
matrice 2D et une matrice 3D (10% du temps de l’entier nieurs système laissent généralement le choix de la version à
du code est utilisé pour cette opération). l’utilisateur. C’est ainsi que sur la machine Pleiades [4], pas
D’un point de vue plus précis, il est clair que l’exécution moins de 5 versions du compilateur Intel, 3 versions de la
du code sur quatre cœurs d’un nœud Woodcrest est environ bilbiothèque MKL sans compter le compilateur libre gfortran
20% plus rapide que l’auto-parallélisation, même optimisée sont disponibles pour les utilisateurs.
manuellement. Cela signifie qu’après avoir aidé le compila-
teur (déroulement des boucles, réordonnancement d’index,
etc.), il est possible d’augmenter la performance d’au moins Conclusion
20%. Afin d’être complet, il a été dit plus haut que la version
parallèle du code prenait 11.07 [s/iter] compilé avec la version On a vu que parvenir à une performance aussi proche
9.1e et 1.99 [s/iter] avec la version 9.0; ceci sur 4 cœurs. que possible de la performance maximum d’un nœud multi-
Cela signifie que cette version du code tourne sur 4 cœurs processeurs multi-cœurs demande du travail. Il est illusoire de
aussi rapidement que la version optimisée manuellement se dire que le compilateur y arrivera seul. Malgré les efforts
sur un seul cœur ! constants des compagnies qui développent des compilateurs
Avant de conclure, un dernier aspect doit être abordé. de plus en plus performants, tout utilisateur/développeur
Afin de se rapprocher le plus possible de la performance HPC devrait être capable de fournir un code source à la per-
maximale autorisée sur un cœur, nous avons observé que formance identique quelle que soit l’option de compilation,
certaines routines BLAS 1 hautement optimisées par Intel en d’autres termes, un code source optimisé à la main pour
(la bibliothèque MKL) ne supportaient pas la scalabilité ! un type d’architecture.

FI 3 – 27 mars 2007 – page 11

Optimisation et parallélisme avec OpenMP: l'approche grain fin

En résumé, la méthodologie présentée dans cet article Références

se décompose en:
I Optimisation sur un processeur [1] Histoire des micro-processeurs, site Web, 2007, www.
(a) Compilation avec l’optimisation la plus basse cpu-info.com
(b) Modification de l’ordre des index [2] Liste des 500 superordinateurs les plus puissants de la
(c) Déroulement des boucles planète, site Web, 2007, www.top500.org
(d) Vectorisation [3] Site officiel du développement d’OpenMP, site Web,
II Parallélisation 2007, www.openmp.org
(a) OpenMP [4] Site du projet Pleiades, site Web, 2007, pleiades.epfl.
(b) ev. Optimisation maximum ch
(c) ev. auto-parallélisation [5] Performance API, site Web , 2007, icl.cs.utk.edu/papi/
Et que ce même utilisateur/développeur garde toujours index.html
en mémoire que la magie n’existe pas; qu’un code médiocre [6] Gruber, R., Volgers, P., De Vita, A., Stengel, M., Tran,
mono-processeur ne donnera jamais de bonnes performances T.-M. , Parameterisation to tailor commodity clusters to
sur plusieurs processeurs. applications. Future Generation Computer Systems,
19:111–120, 2003.
[7] Anderson D.V., Cooper W.A., Fu G.Y., Gengler M.,
Gruber R., Merazzi S., Schwenn U., TERPSICHORE. A
Three-Dimensional Ideal Magnetohydrodynamic Stability
Program, EPFL - Supercomputing Review 3, 29 - 32 ,
1991
[8] W. A. Cooper, S. P. Hirshman, S. Merazzi and R, Gruber,
3D Magnetohydrodynamic Equilibria With Aniostro-
pic Pressure, Computer Physics Communications 72
(1992), 1-13.
[9] Intel Fortran Compiler: Optimizing Applications, Docu-
ment Number 307781-003US, Intel Fortran Compiler
version 9.1 Documentation, Online ftp://download.
intel.com/support/performancetools/fortran/linux/v9/
optaps–for.pdf
On ne fait pas une bonne soupe avec des légumes avariés [10] Intel VTune, Available for download on Intel site, www.
intel.com
[11] OpenMP Application Program Interface Specifications,
Machines disponibles à l’école Version 2.5, May 2005, Available online: www.openmp.
org/drupal/mp-documents/spec25.pdf
Le tab. 11 décrit les différentes machines SMP ou NUMA [12] LAPACK library, OnLine on www.netlib.org
disponibles à l’école. Leur utilisation passe par une politi- [13] Etienne Gondet, Pierre-François Lavallée, Cours
que qui varie selon la machine. Les auteurs conseillent aux OpenMP – CNRS-IDRIS, Version 1.3, 2000, Online,
utilisateurs de se renseigner auprès des ingénieurs système www.idris.fr/data/publications/openmp.pdf
attitrés aux machines. [14] Haldenwang,P. and Labrosse, G. and Abboudi, S.A.
Le tab. 11 présente les machines sur lesquelles il est and Deville, M. , Chebyshev, 3D Spectral and 2D Pseu-
possible d’exécuter un code implémenté avec les directives dospectral Solvers for the Helmholtz Equation, Journal of
de compilation OpenMP. Notons que la performance d’une Computational Physics , Vol 55, pp. 115-128 , 1984
machine NUMA (Accès non uniforme de la mémoire) peut [15] Leriche, E. , Direct Numerical Simulation of a Lid-Dri-
varier en fonction de l’utilisation de cette dernière. ven Cavity Flow by a Chebyshev Spectral Method, PhD
Thesis Numéro 1932, École Polytechnique Fédérale de
Lausanne, 1999 n
Machine Type #Nœuds #cœurs #threads max
per node
Pleiades2+ SMP 99 4 4
Mizar Cluster SMP 224 2 2
Alcor SMP 24 4 4
Mizar ICT NUMA 8 16 -
Altix MX NUMA 1 16 -
Tab. 11 – Machines connues disponibles à l’école

FI 3 – 27 mars 2007 – page 12