Benoit Semelin 2005

Méthodologie M2

Programmation parallèle
pour
le calcul scientifique
Benoît Semelin

2014
Benoit Semelin 2005

Préliminaire
“1000 x ε ~ ε”

- Loguez vous sur [Link] puis sur [Link]

- Ecrire un programme qui:

- Déclare 2 tableaux de doubles tridimensionnels: 512x512x512
- En remplit un: array1(i,j,k)=double(i+j+k)
- Et remplit l’autre ainsi:
i+2 j+2 k+2
array2(i, j, k) = ∑ ∑ ∑ array1(i', j ', k ')
i'=i−2 j '= j−2 k '=k−2

- On se limitera à 2<i<511, 2<j<511, 2<k<511

- Afficher array2(300,300,300)

- Exécuter le programme avec la commande « time ./mon-exécutable »

Benoit Semelin 2005
Benoit Semelin 2005

Introduction

Généralités
Benoit Semelin 2005

Documents en ligne
Ce cours:
[Link]
OpenMP:
[Link] : spécifications officielles

[Link]
[Link]
[Link] : Cours de l'IDRIS

MPI:
[Link] Cours et memo de l'IDRIS
[Link] Standards mpi
Benoit Semelin 2005

Le calcul parallèle:
qu'est ce que c'est?

« Faire coopérer plusieurs processeurs pour réaliser un calcul »

Avantages:
Rapidité:
Pour N processeurs, temps de calcul divisé par N, en théorie…

Taille mémoire:
Pour N processeurs, on dispose de N fois plus de mémoire (en général)

Difficultés:
Il faut gérer le partage des tâches.

Il faut gérer l'échange d'information. (tâches non-indépendantes)

Benoit Semelin 2005

Qu'est ce qui n'est pas

du calcul parallèle.

Quand les processeurs ne coopèrent pas:

Calculs monoprocesseurs séquentiels

Calculs multiprocesseurs:
ü Exécution pour une série de conditions initiales différentes.
ü Problème divisible en sous-problèmes indépendants:
Exemple: Mouvement de N particules test dans un champ
extérieur.

Une architecture parallèle efficace coûte cher, il faut l'utiliser à bon

escient.
Benoit Semelin 2005

Modèles de parallélisme
Architecture matérielle: SISD SIMD MIMD
Single Single Instruction Multiple Instruction
Instruction Multiple Data Multiple Data
Single Data
Architecture Architecture
PC Vectorielle, parallèle
monoprocesseur MMX,SSE, GPU
multiprocesseur
Modèle de programmation:

SPMD MPMD
Le plus utilisé Single Program Multiple Program
Multiple Data Multiple Data
Un seul programme Création dynamique de
(n° processeur = variable) process
Modèle maître-esclave.
Outils de parallélisation:
OpenMP MPI
ordinateur à ordinateur à
mémoire partagée mémoire distribuée
Benoit Semelin 2005

Mémoire partagée / distribuée

Messages

CPU CPU CPU CPU CPU CPU

Écriture

Écriture
Lecture

Lecture
Mémoire Mémoire Mémoire
Mémoire RAM RAM RAM RAM

Mémoire partagée (SMP) Mémoire distribuée

Tous les processeurs ont accès à Chaque processeur possède sa propre

mémoire. Il n'a pas accès à celle des autres.
l'ensemble de la mémoire.
èIl faut gérer l'échange de messages (surcoût)
èAttention aux conflits. èArchitecture bon-marché.

èTrès peu de surcoût de parallélisation. èIl faut ajouter un réseau de com performant.

èLe plus souvent nb proc < 64. èNb de proc ~ illimité.

èArchitecture coûteuse.
Benoit Semelin 2005

Revue d'effectif
dans le monde
Tous les ans, une liste des 500 plus gros ordinateur est publiée sur: [Link]

Rank Name Computer Manufacturer Country Total Cores Rmax (Pflops)

Intel Xeon E5-2692 12C
Tianhe-2
1 2.200GHz, TH Express-2, NUDT China 3120000 33,86
(MilkyWay-2)
Intel Xeon Phi 31S1P
Cray XK7 , Opteron 6274
2 Titan 16C 2.200GHz, NVIDIA Cray Inc. United States 560640 17,59
K20x
BlueGene/Q, Power BQC
3 Sequoia IBM United States 1572864 17,17
16C 1.60 GHz, Custom
K computer, SPARC64
4 K computer Fujitsu Japan 705024 10,51
VIIIfx 2.0GHz
BlueGene/Q, Power BQC
5 Mira IBM United States 786432 8,58
16C 1.60GHz, Custom
Cray XC30, Xeon E5-2670
6 Piz Daint 8C 2.600GHz, , NVIDIA Cray Inc. Switzerland 115984 6,27
K20x
Xeon E5-2680 8C
7 Stampede 2.700GHz, Infiniband FDR, Dell United States 462462 5,16
Intel Xeon Phi SE10P
BlueGene/Q, Power BQC
8 JUQUEEN IBM Germany 458752 5,00
16C 1.600GHz
BlueGene/Q, Power BQC
9 Vulcan IBM United States 393216 4,29
16C 1.600GHz

10 CS-storm Cray, Intel Xeon E5-2660v2 IBM United States 72800 3.57
bullx DLC, Xeon E5-2690v3
26 Occigen Bull SA France 50544 2.1
12C 2.6GHz, Infiniband FDR
Benoit Semelin 2005

Revue d'effectif
pour la recherche en France

GENCI: Grand Equipement National de Calcul Intensif

- le CCRT: Tera 100 (Curie): 1.35 Pétaflops

- le CINES: Bull (Occigen) 2.1 Pflops

- l’ IDRIS: IBM Bluegene (Turing) 1.25 Pflops + IBM (Ada) 230 Tflops

Une procédure de demande d’heures de calcul unique, sur projet

détaillé.

PRACE: - Demande à l’échelle européene.

- ~50. 106 h de calcul par demande.
- Le code doit « bien » fonctionner sur 4000 cœurs.
Benoit Semelin 2005

Revue d'effectif 3:
à l’obsevatoire
Dans le cadre de cette option, nous travaillerons sur:

MesoPSL:

1472 cœurs, 4Go/cœurs, 64 Go par nœud.

Accès sur projet (demande ~ 1 page actuellement)

Momentum:

Acutellement: 1400 cœurs AMD, 64 cœurs/nœud, 128 Go/nœud, 8 GPU Fermi.

Accès fermé. Cluster réservé à un projet.

Evolution
Benoit Semelin 2005

Outils de parallélisation
La parallélisation peut-être effectuée à divers niveaux:
Langages, ou extensions de langages:
è CUDA
è OpenCL
Bibliothèques:
è Message Passing Interface (MPI)
è Pthreads (langage C)
Directives de compilation:
è OpenMP
è Directives d’accélération pour GPU ou autre.
Compilateurs: efficacité très faible.
è Intel Fortran/C compiler: gratuit.
Benoit Semelin 2005

OpenMP
Benoit Semelin 2005

OpenMP:

modèle de programmation
Benoit Semelin 2005

La parallélisation facile:
OpenMP
OpenMP est un ensemble de directives de compilation pour
paralléliser un code sur une architecture SMP (interfaces Fortran, C
et C++)
Le compilateur interprète les directives OpenMP (si il en est capable!)
Les standards d'OpenMP datent de 1997, ceux d'OpenMP-2 de 2000, OpenMP-3
2008. Les développeurs de compilateurs les implémentent.
Modèles d'exécution OpenMP: Thread ≃ process (≠
Procces processeur)
princip
al

Exécution Exécution Exécution Exécution Exécution

séquentielle parallèle séquentielle parallèle séquentielle
monoprocess à 9 process monoprocess à 9 process monoprocess
Benoit Semelin 2005

Comment « faire tourner »

un code OpenMP
Voici quelques points de repère:

Écrire le code séquentiel (tout ou partie)

Débuguer le mieux possible
Se placer sur une machine SMP multiprocesseur (optionnel)
Insérer les directives OpenMP
Compiler avec l'option appropriée: ifort -openmp toto.f90 -o [Link]
Définir le nombre de CPU: export OMP_NUM_THREADS =4
Lancer normalement: ./[Link]
Débuguer...
Évaluer les performances, en particulier le « speedup »:
Temps exécution 1 procs / Temps exécution N proc
Benoit Semelin 2005

Syntaxe d'une directive

OpenMP

Voici un exemple de directive pour f90:

Espace en 6ième colonne pour f77

!$OMP_DO_SCHEDULE(DYNAMIC,500)

Sentinelle: Directive: Clause:

Doit être le premier Une directive par ligne. Optionnelle. Modifie le
caractère non-blanc de la comportement de la
ligne. « !$ » peut servir de directive. Il peut y en
sentinelle de compilation avoir plusieurs, elle peut
conditionnelle pour une être répétée, l'ordre est
ligne de code normale. indifférent

Syntaxe en C/C++: #pragma_omp_for_schedule(dynamic,500)

Benoit Semelin 2005

OpenMP:

régions parallèles,
comportements des variables
Benoit Semelin 2005

Définir une région parallèle

SUBROUTINE compute_grad(field,grad,n,size,TYPEMIN)
Comment amorcer une zone
USE VARIABLES
d'exécution multiprocesseur:
INTEGER, INTENT(IN) :: n
REAL(KIND=8), INTENT(IN) :: size
REAL(KIND=8), DIMENSION(n), INTENT(IN) :: field
Ø Début de zone: REAL(KIND=8), DIMENSION(n), INTENT(OUT) :: grad
!$OMP PARALLEL INTEGER, INTENT(IN) :: typemin
REAL(KIND=8) :: val,g1,g2,fx
#pragma omp parallel {code} INTEGER :: i,ip,im

Ø Fin de zone: !$OMP PARALLEL

! Autre directive OMP
!$OMP END PARALLEL do i=1,NCELLNOW
Dans la zone correspondante, N threads ip=i+1
im=i-1
sont exécutés en parallèle. N est fixé par if(ip == NCELLNOW+1) ip=1
if(im == 0) im=ncellnow
la variable d'environnement if(celltype(i) >= typemin) then
g1=(field(i)-field(im))/(cellsize(i)+cellsize(im))*2.
OMP_NUM_THREADS. g2=(field(ip)-field(i))/(cellsize(i)+cellsize(ip))*2.
fx=(cellsize(im)+cellsize(i))/(cellsize(im)+2.*cellsize(i)+cellsize(ip))
grad(i)=(1.-fx)*g1+fx*g2
Que devient la valeur d'une variable à endif
enddo
l'entrée et à la sortie d'une zone !$OMP END PARALLEL
parallèle?? Que vaut-elle sur 2
END SUBROUTINE compute_grad
processeurs différents ?
Benoit Semelin 2005

Comportements possibles
d'une variable
Variable SHARED Variable PRIVATE

X=1. X=1.

X=1. X=1. X=1. X

X=1.
1=? X
X=1.
2=? X
X=1.
3=? Xi indéterminé.
... ... ... ... ... ... contrôlable par
... ... ... ... ... ... FIRSTPRIVATE
... X=X+1. ... ... X=X+1.
... ...
X=2. Délai ... X=2.
... ... ...
(cache, registre) X=2. X1=2.53 X2=124. X
X=2.
3=2.
contrôlé par
FLUSH

X=2. X=2.
X=? X indéterminé.
contrôlable par
LASTPRIVATE
Benoit Semelin 2005

Comportement SHARED/PRIVATE
d'une variable
SUBROUTINE compute_grad(field,grad,n,cellsize,TYPEMIN)
SHARED et PRIVATE sont deux
USE VARIABLES
clauses qui spécifient le
INTEGER, INTENT(IN) :: n
comportement des variables dans REAL(KIND=8), INTENT(IN) :: cellsize
REAL(KIND=8), DIMENSION(n), INTENT(IN) :: field
une zone parallèle: REAL(KIND=8), DIMENSION(n), INTENT(OUT) :: grad
INTEGER, INTENT(IN) :: typemin
Ø SHARED: la variable a la même REAL(KIND=8) :: val,g1,g2,fx
INTEGER :: i,ip,im
valeur sur tout les processeurs. C'est le
statut par défaut. Attention à la !$OMP PARALLEL PRIVATE(ip,im,g1,g2,fx)
synchro si on la modifie sur un ! Autre directive OMP

processeur (directive FLUSH). do i=1,NCELLNOW

ip=i+1
im=i-1
Ø PRIVATE: la variable a une valeur if(ip == NCELLNOW+1) ip=1
if(im == 0) im=ncellnow
différente sur chaque processeur if(celltype(i) >= typemin) then
(espace mémoire correspondant g1=(field(i)-field(im))/(cellsize(i)+cellsize(im))*2.
g2=(field(ip)-field(i))/(cellsize(i)+cellsize(ip))*2.
dupliqué) et indéterminée en entrée de fx=(cellsize(im)+cellsize(i))/(cellsize(im)+2.*cellsize(i)+cellsize(ip))
zone. C'est le statut de toute variable grad(i)=(1.-fx)*g1+fx*g2
endif
déclarée à l'intérieur d'une zone enddo
parallèle. !$OMP END PARALLEL

END SUBROUTINE compute_grad

Benoit Semelin 2005

Autres clauses de
comportements de variables
On peut changer le comportement défaut:
!$OMP PARALLEL DEFAULT(PRIVATE | SHARED | NONE)

Si on utilise NONE, on doit préciser le comportement de toutes les variables qui

apparaissent dans la zone parallèle.

Clauses supplémentaires:
Ø FIRSPRIVATE(X): X est PRIVATE et initialisée à sa valeur juste avant la zone
parallèle.
Ø REDUCTION(*,X): X est PRIVATE et un produit des valeurs de X sur les
différents threads est effectué en fin de zone parallèle et stocké dans X. L'opérateur
peut-être + , - , * , .OR. , .AND. , MAX, MIN, etc... Remplace les statuts PRIVATE
ou SHARED.
Ø LASTPRIVATE(X): X est PRIVATE. La valeur de X du thread exécutant la
dernière mise à jour de X est conservée en fin de zone parallèle.
Ø COPYPRIVATE(X): Pour diffuser une variable privée (Directive SINGLE
uniquement).
Benoit Semelin 2005

Variable PRIVATE mais

globale... dans la zone paralèlle
Ø Une variable PRIVATE est locale.
Ø Comment conserver la valeur d'une variable PRIVATE d'une zone parallèle à
l'autre ?
Ø Il est intéressant de pouvoir appeler une ou plusieurs procédure/fonction dans une

zone parallèle.
Ø Certaines procédures/fonctions utilisent et modifient des variables globales.

!$OMP THREADPRIVATE(X,/COORD/) (Dans le var.f90)

La variable X et le bloc commun /COORD/ seront PRIVATE mais globaux dans

chaque thead des régions parallèles.

!$OMP PARALLEL COPYIN(/COORD/) (Dans le .f90)

Les valeurs de /COORD/ sont copiées dans les répliques privées de chaque thread en
entrée dans les régions parallèles.
Benoit Semelin 2005

OpenMP:

partage du travail
entre threads
Benoit Semelin 2005

Partage du travail:
distribuer une boucle
SUBROUTINE compute_grad(field,grad,n,cellsize,TYPEMIN)
La directive DO se place juste avant le
USE VARIABLES
début d'une boucle, elle répartit les
INTEGER, INTENT(IN) :: n
itérations entre les processeurs. REAL(KIND=8), INTENT(IN) :: cellsize
REAL(KIND=8), DIMENSION(n), INTENT(IN) :: field
REAL(KIND=8), DIMENSION(n), INTENT(OUT) :: grad
è Pas de DO WHILE ! INTEGER, INTENT(IN) :: typemin
REAL(KIND=8) :: val,g1,g2,fx
è Le mode de répartition dépend de la INTEGER :: i,ip,im

clause optionnelle SCHEDULE !$OMP PARALLEL PRIVATE(ip,im,g1,g2,fx)

!$OMP DO SCHEDULE(DYNAMIC,20)
è Le mode de répartition par défaut do i=1,NCELLNOW
ip=i+1
dépend de l'implémentation im=i-1
if(ip == NCELLNOW+1) ip=1
d'OpenMP. if(im == 0) im=ncellnow
if(celltype(i) >= typemin) then
Les itérations doivent être g1=(field(i)-field(im))/(cellsize(i)+cellsize(im))*2.
g2=(field(ip)-field(i))/(cellsize(i)+cellsize(ip))*2.
indépendantes (ordre indifférent). fx=(cellsize(im)+cellsize(i))/(cellsize(im)+2.*cellsize(i)+cellsize(ip))
grad(i)=(1.-fx)*g1+fx*g2
A la fin de la boucle on insère la endif
enddo
directive END DO. Les threads !$OMP END DO
!$OMP END PARALLEL
se synchronisent (on attend la
END SUBROUTINE compute_grad
dernière).
Benoit Semelin 2005

Stratégies de répartition
des itérations
La clause SCHEDULE permet de contrôler la stratégie de répartition des itérations
d'une boucle entre les processeurs, elle admet 2 arguments: SCHEDULE(stratégie,N).

Stratégies possibles:
ü STATIC: Chaque thread reçoit à tour de rôle N itérations à traiter. La distribution
s'effectue dans un ordre fixé, éventuellement en plusieurs tours. N est optionnel, si il n'est
pas spécifié N~ nb itération / nb threads.
ü DYNAMIC: Chaque thread reçoit N itérations à traiter. Dès qu'un thread a fini, il en
reçoit N autres, jusqu'à épuisement du travail. N est optionnel, par défaut N=1.
ü GUIDED: Les itérations sont divisées en paquets de taille exponentiellement
décroissante. Les paquets sont distribués dynamiquement. La taille du plus petit paquet
est N.
ü RUNTIME: Le choix de la stratégie est reporté au moment de l'exécution. Il sera
alors déterminé par le contenu de la variable d'environnement OMP_SCHEDULE.

DYNAMIC et GUIDED assurent un meilleur équilibrage de charge que STATIC.

Benoit Semelin 2005

Un exemple:
calcul de potentiel périodique
PROGRAM POTPER
IMPLICIT NONE
Programme pour calculer une REAL(KIND=8),PARAMETER :: BOXSIZE=1.
INTEGER, PARAMETER :: NREP=200
interaction gravitationnelle avec INTEGER, PARAMETER :: NDEF=40
INTEGER :: i,j,k,l
des conditions de bord REAL(KIND=8) :: x1,y1,z1,x2,y2,z2,pot,dist

périodiques. y1=0.
z1=0.
do i=1,ndef
x1=dble(i)/ndef*BOXSIZE
pot=0.
Remarques: !$OMP PARALLEL PRIVATE(k,l,x2,y2,z2,dist) REDUCTION(+:pot)
!$OMP DO SCHEDULE(DYNAMIC)
Il existe un algorithme plus do j=-nrep,nrep
x2=dble(j)*BOXSIZE
efficace (décomposition de la do k=-nrep,nrep
somme en deux parties). y2=dble(k)*BOXSIZE
do l=-nrep,nrep
On pourrait paralléliser la z2=dble(l)*BOXSIZE
dist=sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2)
première boucle! if(dist < BOXSIZE*NREP*0.95) then
Remarquer l'utilisation de la pot=pot+(x2-x1)/(sqrt((x2-x1)**2+(y2-y1)**2+(z2-z1)**2))**3
endif
clause REDUCTION. enddo
Zones parallèles longues => bon enddo
enddo
speedup !$OMP END DO
!$OMP END PARALLEL
print*,x1,pot
enddo
END PROGRAM POTPER
Benoit Semelin 2005

Partage du travail
en l'absence de boucle
La directive SECTIONS amorce une
zone où le code est découpé en ...
morceaux. END SECTIONS clôt cette
!$OMP PARALLEL
zone !$OMP SECTIONS
Ces morceaux sont séparés par des
directives SECTION. !$OMP SECTION
CALL UPDATE_XPOS()
Chaque morceau (section) sera exécuté
une fois unique par un des threads. !$OMP SECTION
CALL UPDATE_YPOS()
L'ordre d'exécution des sections doit
!$OMP SECTION
être indifférent! CALL UPDATE_ZPOS()

SECTIONS admet PRIVATE, !$OMP END SECTIONS

!$OMP END PARALLEL
FIRSTPRIVATE, LASTPRIVATE, et
...
REDUCTION comme clauses.
Benoit Semelin 2005

Partage du travail avec

WORKSHARE
WORKSHARE, une directive trop
générale/ambitieuse? !$OMP PARALLEL
!$OMP WORKSHARE

D'après les spécifications d'OpenMP X(1:N) = sin(THETA(1:N))

WORKSHARE doit diviser le travail de telle Y(1:N) = cos(THETA(1:N))
manière que chaque instruction de la zone soit WHERE (Y .ne. 0) P = X / Y
exécutée exactement 1 fois dans un thread, en
respectant la sémantique du code...Si les
FORALL( i = 1 : N , j = 1 : N , i /= j )
instructions ne sont pas indépendantes, pas de
pot(i)=pot(i)+1./(x(i)-x(j))
speedup.
END FORALL
Aux fabricants de compilateurs de l'implémenter.
!$OMP END WORKSHARE
!$OMP END PARALLEL
En pratique, utilisable en f90 pour:
Opérations sur les tableaux.
Instruction WHERE
Instruction FORALL
WORKSHARE n'admet AUCUNE clause.
Benoit Semelin 2005

Exécution exclusive:
directives SINGLE et MASTER
Si, dans une zone parallèle, on souhaite qu'une partie du code soit
exécutée par un seul thread, on a le choix entre 2 directives:

!$OMP SINGLE / !$OMP END SINGLE :

Le premier thread qui aborde la zone l'exécute. Les autres sautent à la fin et
attendent que celui qui exécute ait fini. SINGLE admet les clauses
PRIVATE et FIRSTPRIVATE. END SINGLE admet NOWAIT et
COPYPRIVATE.
!$OMP MASTER / !$OMP END MASTER :
Le thread « master » (numéro 1) exécute la zone. Les autres sautent la zone,
et n'attendent pas! MASTER n'admet aucune clause.

Ces directives peuvent par exemple servir à faire des entrées/

sorties.
Benoit Semelin 2005

OpenMP:

synchronisation des threads

Benoit Semelin 2005

Problèmes de
synchronisation.

Dans une zone parallèle, l'utilisation de variables SHARED peut

créer des problèmes de synchronisation. Exemples:

3 threads calculent X(1),X(2) et X(3), puis chacun utilise les 3 valeurs. Il

faut s'assurer que les autres ont finit leur calcul.

N threads incrémentent la même variable. Il faut s'assurer qu'ils ne tentent

pas d'écrire en même temps dans l'emplacement mémoire de la variable.

Problème de vidage de registre et de tampon.

Benoit Semelin 2005

Synchronisation simple:
BARRIER
La directive BARRIER synchronise les threads: tous s'arrêtent au
niveau de la directive jusqu'à ce que le dernier soit arrivé. Puis ils
continuent tous. Syntaxe:
!$OMP BARRIER
N'admet aucune clause.
Une BARRIER est implicitement incluse dans les directives
suivantes:
Ø END PARALLEL
Ø END DO (sauf clause NOWAIT)
Ø END SECTIONS (sauf clause NOWAIT)
Ø END SINGLE (sauf clause NOWAIT)
Ø END WORKSHARE (sauf clause NOWAIT)
Pas de BARRIER implicite de END MASTER !
Benoit Semelin 2005

Mettre à jour la mémoire:

directive FLUSH
!$OMP FLUSH(x,y)
Lorsqu'un thread rencontre une directive FLUSH, il met à jour les
variables spécifiées entre parenthèses, c'est à dire:
ü Il vide les tampons d'écriture, les caches, les registres.

ü Il vérifie si un autre thread n'a pas modifié la variable.

FLUSH n'a d'utilité que pour les variables visibles des autres threads.
Si aucune variable n'est spécifiée, toutes les variable visibles sont mises à jour.

La directive FLUSH est implicitement appelée par les directives suivantes:

BARRIER, CRITICAL, END CRITICAL, END DO , END SECTIONS , END
SINGLE, END WORKSHARE, ORDERED, END ORDERED, PARALLEL et
END PARALLEL.

On a rarement besoin d'utiliser FLUSH directement. Il est plus simple (et plus lent)
d'utiliser BARRIER.
Benoit Semelin 2005

Éviter les conflits:

directives ATOMIC et CRITICAL
Comment éviter que 2 threads tentent de changer en même temps la
même variable shared.
!$OMP PARALLEL PRIVATE(XLOCAL)
Directive ATOMIC: !$OMP DO
do I=1,N
La ligne de code qui suit la directive ATOMIC et qui call compute(xlocal(i),i)
modifie une variable X, est exécutée atomiquement,
!$OMP ATOMIC
c'est à dire jamais simultanément pas deux threads. xglobal = xglobal + xlocal(i)
(Voir référence OpenMP pour la forme de la ligne de
enddo
code concernée). !$OMP END DO
!$OMP END PARALLEL
ATOMIC peut être plus performant que CRITICAL.
!$OMP PARALLEL PRIVATE(XLOC,YLOC)

Directive CRITICAL: call compute_xloc()

Même fonction que ATOMIC, mais concerne une call compute_yloc()

zone de code. Un seul thread peut accéder à cette !$OMP CRITICAL(ZONE1)

xglob=xglob+xloc
zone simultanément. La zone se termine par END yglob=yglob+yloc
CRITICAL. Si on a plusieurs zones CRITICAL, il !$OMP END CRITICAL(ZONE1)
faut les nommer pour les rendre indépendantes. !$OMP END PARALLEL
Benoit Semelin 2005

Dans l'ordre:
directive ORDERED
La directive ORDERED permet, à !$OMP PARALLEL PRIVATE(XLOCAL)

l'intérieur d'une boucle parallélisée, !$OMP DO ORDERED SCHEDULE(DYNAMIC)

d'exécuter une zone séquentiellement, do I=1,N
call compute(xlocal(i),i)
c'est à dire thread par thread, dans !$OMP ORDERED
write(*,*) i, xlocal(i)
l'ordre des indices croissant. !$OMP END ORDERED
enddo
- Cela permet de faire une entrée-sortie !$OMP END DO

ordonnée dans une zone parallèle. !$OMP END PARALLEL

- C'est un outil de débugage. Cela permet de

vérifier l'indépendance des itérations d'une
boucle.

- Si la zone ORDERED représente une

fraction du temps de calcul d'une itération,
supérieure à 1./OMP_NUM_THREAD, cela
ralentit l'exécution.
Benoit Semelin 2005

OpenMP:

variables d'environnement
et bibliothèque standard
Benoit Semelin 2005

Variables d'environnement

Il existe 4 variables qui définissent l'environnement OpenMP:

u OMP_DYNAMIC: booléen. Permet l'ajustement dynamique du nombre de threads.
Si TRUE, l'utilisateur déclare le nombre maximal de thread, et le système en donne
un nombre inférieur ou égal, en fonction des disponibilités.
u OMP_NUM_THREADS: entier. Fixe le nombre (maximal) de threads.
u OMP_NESTED: booléen. Permet d'imbriquer les régions parallèles les unes dans
les autres: chaque threads peut se subdiviser. Délicat à utiliser. Inutile pour le calcul
scientifique?
u OMP_SCHEDULE: chaîne de caractères. Spécifie la stratégie de distribution des
itérations de la directive DO associée à la clause RUNTIME.
Exemples:
export OMP_DYNAMIC=TRUE
export OMP_NUM_THREADS=8 (souvent indispensable)
export OMP_NESTED=FALSE
export OMP_SCHEDULE=''GUIDED 4 ''
Les valeurs des variables peuvent être modifiées pendant l'exécution par des fonctions
de la bibliothèque standard.
Benoit Semelin 2005

Bibliothèque standard:
contrôle de l'environnement
On peut souvent se passer complètement de ces fonctions.
Modèles d'exécution: (appel depuis une zone séquentielle)
è OMP_SET_DYNAMIC(boolean): Subroutine. Active/désactive l'ajustement
dynamique du nombre de thread: fixe OMP_DYNAMIC.
è OMP_GET_DYNAMIC(): Function. Retourne la valeur actuelle de OMP_DYNAMIC.
è OMP_SET_NESTED(boolean): Subroutine. Fixe OMP_NESTED
è OMP_GET_NESTED(): Function. Retourne la valeur actuelle de OMP_NESTED .

Contrôle du nb de thread/processeur:
è OMP_SET_NUM_THREADS(entier): Subroutine. Fixe le nombre de threads
(maximum) pour les prochaines zones parallèles (valeur de OMP_NUM_THREADS).
Appel depuis une zone séquentielle.
è OMP_GET_NUM_THREADS(): Function. Retourne le nb réel de threads utilisées à
l'instant t.
è OMP_GET_MAX_THREADS(): Function. Retourne la valeur de
OMP_NUM_THREADS (nb maximal de threads).
è OMP_GET_NUM_PROCS(): Function. Retourne le nombre de processeurs utilisés.
Benoit Semelin 2005

Bibliothèque standard:
contrôle manuel de la parallélisation

Exécution conditionnelle/ partage des tâches manuel:

è OMP_IN_PARALLEL(): Boolean function. Détermine si on est dans une région
parallèle. La directive PARALLEL peut comporter une clause de condition (cf
spécifications OpenMP). Utile dans ce cas.
è OMP_GET_THREAD_NUM(): Fonction entière. Retourne le numéro du thread.
Permet, par exemple, de faire le partage du travail ''à la main'', sans utiliser de directive
OpenMP ( Masochiste sur une architecture SMP).

Utilisation de verrous:
Un verrou est libre ou possédé par un thread. Une série de fonction OpenMP
permet de les manipuler, par exemple d'attendre à un point du code que le
verrou soit libéré par un autre thread.

Permet des comportements similaires à la directive CRITICAL, mais

avec un contrôle plus fin. Voir spécifications OpenMP.
Benoit Semelin 2005

Parallélisation sur architectures

à mémoire distribuée:

décomposition en domaines
Benoit Semelin 2005

Décomposition en
domaine / fonction
Schéma basique d'un ordinateur: x, y
CPU
z Mémoire
- Une unité de calcul opère sur des données z=x+y
qu’elle lit et écrit en mémoire.

Si on a N CPU et N mémoires, il faut:

Ø Répartir les calculs:
Chaque unité de calcul va assurer un type de calcul différent sur l'ensemble des
données. C'est le principe de décomposition par fonction.
Ø Répartir les données:
Toutes les unités de calcul effectuent les mêmes calculs sur la partie de données qui
leur est attribuée. L’accès aux données des autres unités de calcul se fait par échange
de messages. C'est le principe de la décomposition en domaines.

La séparation est arbitraire, en général on combine les deux. Mais,

pour le calcul scientifique, c'est la décomposition en domaines qui est
cruciale.
Benoit Semelin 2005

Décomposition en domaine:
stratégie.
Les contraintes:
- Les communications (échange de données entre domaines) ont un coût en temps
CPU: il faut minimiser la quantité de communication (nombre et taille).
- Il faut que chaque domaine représente la même quantité de travail. Il faut minimiser
les fluctuations de charge et assurer l'équilibrage dynamique des charges.

Une décomposition idéale:

Aucun besoin de communication entre domaines.
Domaines statiques: charge constante, équilibré.

Ce n'est pas du parallélisme!

Une décomposition satisfaisante:
Temps de communication << temps de calcul.
Equilibrage avec variation relative de charge < 10%.
Temps de redéfinition dynamique des domaines << temps de calcul.
Benoit Semelin 2005

Décomposabilité

La possibilité/efficacité d'une décomposition en domaines dépend du

problème-algorithme. Exemples:

Pb-algo non décomposable (non parallélisable...):

Intégration d'une ODE par méthode Euler, Runge-Kutta...
Dynamique un système à trois corps (cf ODE).

Pb-algo à décomposition délicate ou inefficace:

Dynamique d'un système N corps.
Résolution d'une EPD sur grille adaptative et pas de temps adaptatif (équilibage).
Réseau de réactions chimiques (vitesses de réaction).

Pb-algo à décomposition simple et efficace:

Résolution d'une EDP sur grille non adaptative.
Simulation Monte-Carlo pour certains pbs (équilibrage de charge simple)
Benoit Semelin 2005

Exemple de décomposition 1:
résolution d'EDP sur grille.
Résolution de l'équation de Burger: ∂v ∂v 1 ∂ 2v
+v =
∂t ∂x Re ∂x 2
On résout l'équation sur une grille 1-D, par une méthode de différences finies, et avec
un schéma explicite (!) à pas de temps constant.

Taille de cellule constante: Grille adaptative :

domaines statiques domaines dynamiques.

Le domaine 2 se
restreint spatialement
mais conserve autant
de cellules que 1 et 3.
Benoit Semelin 2005

Exemple de décomposition 2:
L'algorithme Treecode.
Dynamique d'un système N-
corps:

Calcul des forces en N ln(N), par

utilisation du développement
multipolaire. Contrôle de la précision θ
par θ, l'angle d'ouverture. C'est un
algorithme dit « en arbre ».

Conséquences pour la parallélisation:

Moins d'interactions entre zones

distantes.

=> Des domaines ''compacts'' dans

l'espace (x,y,z) minimisent les
communications.
Benoit Semelin 2005

Algorithme Treecode:
exemple de décomposition.

Méthode par dichotomie. Décomposition de l'arbre. Tesselation Voronoï.

Simple et efficace. Pour les algo oct-tree. Communications mimales.
Gadget 1 et 2(code public) Construction des domaines
coûteuse.

Plusieurs solutions possibles. Ca dépend des détails de

l'implémentation de l'algorithme commun.
Benoit Semelin 2005

Parallélisation sur architectures

à mémoire distribuée:

communications.
Benoit Semelin 2005

Temps de transmission
d'un message.
Sur une architecture à mémoire distribuée, les unités de calculs doivent
échanger des messages. Ces communications ne sont pas instantanées.

Le coût dépend du réseau de communiations et de l'implémentation

MPI. des contributions au
Détails
coût de communication. voir
graphique.

Conséquences:

Il vaut mieux envoyer 1

gros message que plein de
petits. Il faut grouper les
communication!
Benoit Semelin 2005

Modes de communications
entre process.
Exemple de mécanisme de communication: (two-sided, buffered,
synchrone)

Lecture Ecriture
CPU 1 CPU 1
RAM

RAM
Ecriture ''Prêt à envoyer'' Lecture
Buffer de Buffer ''Prêt à recevoir'' Buffer de Buffer
réception d'envoi réception d'envoi

Transmission de
données
Mode de commnications: définitions
One-sided / Two-sided:
Un seul processeur gère tout les éléments de la com / les deux processeurs intervienent, ils
gèrent chacun une partie des éléments de la com.
Synchrone / Asynchrone:
Le processeur qui envoit (1) attend / n'attend pas que le proccesseur qui recoit (2) ait
signalé qu'il était prêt. En mode asynchrone, le processeur (1) transmet les données qui
sont écrites dans le buffer de réception du processeur (2). Il viendra les lire plus tard,
quand le code lui dira de faire.
Benoît Semelin 2005

Communication globales:
les contraintes matérielles
Dans un ordinateur multi-processeur, chaque processeur n'a pas un lien avec tout les
autres pour les communications (trop cher). Il faut pourtant tenter de l'optimiser.

Deux critères de performance pour un réseau de n processeurs:

- Nombre de liens moyen pour joindre 2 processeurs nl
- ''Bande passante de bissection'' Bb (bande passante entre 2 moitiés du réseau)

Ligne: nL~n/3 , Anneau: nL~n/4 ,

Bb=1 Bb=2
Grille: nL~√n , Bb=√n Tore: nL~√n/2 , Bb=2√n

Hypercube dim D: nL~D , Arbre hiérarchique: nL~ log(n) , Bb=n/

Bb=n/2 2...?
Benoit Semelin 2005

Communications globales
sur l'hypercube
Une implémentation naïve de la réduction d'une variable (somme des valeurs
sur les N processeurs) nécessite N communications successives:
X1+X2+X3 X1+X2 X1
CPU8 CPU7 CPU6 CPU5 CPU4 CPU3 CPU2 CPU1
X8 X7 X6 X5 X4 X3 X2 X1

Une implémentation de la
CPU7 CPU8
réduction en utilisant le principe
1 X7 1 X8
de l'hypercube nécessite log(N)
communications successives. CPU4 CPU3
Le principe de l'hypercube peut servir pour X4 X3
organiser des communications globales
type REDUCTION dans un code. MPI s'en 2 CPU5 CPU6
2
sert pour ses fonctions de communications 1 X5 1 X6
globales génériques. Si ces routines
génériques ne répondent pas aux besoins de
l 'utilisation, il peut utiliser le principe de CPU1 3 CPU2
l'hypercube lui même. X1 X2
Benoit Semelin 2005

Performances d'un
code parallèle.
Coût des communications:
Si le temps de communications = 10% du temps de calcul, inutile d'utiliser
plus de 5-10 processeurs. Il faut minimiser le temps de communication.
Pistes pour réduire les communications:
- Recouvrement des bords de domaines (Simulation sur grille)
- Duplication des données critiques (Treecode)
Equilibrage de charge:
Si un processeur est 1,3 fois plus lent qu'un autre, on obtient le même temps
de calcul qu'avec 0.7 fois moins de processeurs.
Equilibrage dynamique des charges (''overload diffusion''):
- Applicable, par ex, quand on calcule une évolution temporelle.
- Chaque domaine a une taille Si (t) variable . A chaque pas de temps, on
calcule Ti(t), le temps de calcul sur chaque processeur.
- On adapte Si, par exemple: Si (t+dt)=Si (t)*Tmoy(t)/ Ti(t) (stabilité?)
Benoit Semelin 2005

MPI
Benoit Semelin 2005

Introduction
MPI est une bibliothèque de communication pour le parallélisme sur
architectures à mémoire partagée (Fortran, C, C++).
Les standards MPI:
MPI 1 (1994): Le coeur de la bibliothèque
- Communicateurs
- Communications locales et globale
- Définition de type de variable
- Topologies
MPI 2 (1997): extensions et support C++ et f90
- Communications ''one-sided''
- Entrées-sorties parallèles (non traitées ici)

Implémentations gratuites:
LAM: [Link] (remplacé par OpenMPI)
MPI 1 et 2 (en partie) sur grille de PC ou architectures SMP.
MPICH2: [Link]
MPI 1 et 2 (en partie) sur grille de PC ou architectures SMP.
Benoit Semelin 2005

MPI:

fonctions générales,
environnement
Benoit Semelin 2005

Environnement système
program EX_0
Il faut inclure la bibliothèque
implicit none
MPI. include ''mpif.h'' !(ou USE MPI avant implicit none)
integer :: nb_procs,rang,err
- mpif.h en fortran. call MPI_INIT (err)
- mpi.h en C/C++. call MPI_COMM_SIZE ( MPI_COMM_WORLD ,nb_procs,err)
call MPI_COMM_RANK ( MPI_COMM_WORLD ,rang,err)
print *, 'Je suis le processus' ,rang, 'parmi' ,nb_procs
Exemples de commandes call MPI_FINALIZE (err)
de compilation: end program EX_0

Sur tucana:
- mpif90 EX_0.f90 (ou mpif77, mpicc, mpiCC)
- mpd & (la première fois)
- mpiexec -n 8 EX_0.out
- mpdallexit
Benoit Semelin 2005

Initialisation et sortie
program EX_0
Avec MPI, on définit une seule
zone parallèle dans le implicit none
include ''mpif.h'' !(ou USE MPI avant implicit none)
programme, souvant l'ensemble integer :: nb_procs,rang,err
du code.
call MPI_INIT (err)
call MPI_COMM_SIZE ( MPI_COMM_WORLD ,nb_procs,err)
Les parties du code avant et après call MPI_COMM_RANK ( MPI_COMM_WORLD ,rang,err)
print *, 'Je suis le processus' ,rang, 'parmi' ,nb_procs
la zone parallèle sont locale sur call MPI_FINALIZE (err)
chaque processeur. Pas de
end program EX_0
communications possibles.
Début de zone parallèle Fin de zone parallèle
En fortran: call MPI_INIT(err) En fortran: call MPI_FINALIZE(err)
En C: MPI_Init(); En C: MPI_Finalize();
En C++: MPI::Init(); En C++: MPI::Finalize();
Fonction collective! (appelée par tous les process) Fonction collective! (appelée par tous les process)

Pour une définition précise des fonctions (type, arg optionnels, etc...) voir aide
mémoire IDRIS et manuel de référence MPI en ligne.
Benoit Semelin 2005

Goupes et communicateurs:
objets de base de MPI
Groupes et communicateurs sont des objets MPI opaques: on les
manipule par une ''poignée'' (variable de type INTEGER en fortran)
- Groupe: Ensemble de process.
0 1

2 3

- Intracommunicateur: Struture d'échange de messages entre process

contenant un groupe unique.
0 1

2 3

- Intercommunicateur: Struture d'échange de messages entre process contenant

deux groupes.

1
0 2
1 0 1

3
2 4
3 2 3
Benoit Semelin 200

Manipulation des
groupes et communicateur.
L'appel à MPI_INIT() définit un communicateur par défaut: MPI_COMM_WORLD (constante
entière MPI)
Souvant, MPI_COMM_WORLD est le seul usage des communicateurs/groupes que fait un code.
Mais il est parfois intéressant de créer des sous-groupes/sous-communicateur.
Quelques fonctions pour créer-manipuler-détruire les groupes et communicateur.
MPI_COMM_SIZE(comm,size,err): locale. Renvoie size, le nb de process dans comm.
MPI_COMM_RANK(comm,rank,err): locale. Renvoie rank, le n° du process appelant.

MPI_COMM_GROUP(comm,group,err): locale! Crée un groupe group avec tout les process de

comm.
MPI_GROUP_INCL(group,n,rank_array,subgroup,err): locale. Crée un groupe subgroup, avec
un sélection de n process de group définie par rank_array.
MPI_COMM_CREATE(comm,subgroup,subcomm,err): collective. Crée l'intra-communicateur
subcomm de subgroup,
MPI_GROUP_SIZE(group,size,err): locale. Renvoie size, le nb de process dans group.
MPI_GROUP_RANK(group,rank,err): locale. Renvoie rank, le n° du process appelant.

MPI_GROUP_FREE(group,err): locale. Désalloue le groupe group.

MPI_COMM_FREE(comm,err): collective. Désallooue le communicateur comm.
Benoit Semelin 200

MPI:

communications point à point

Benoit Semelin 200

Structure d'un message MPI

Par l'intermédaire d'un message, un process peut envoyer des données à
un autre.

En plus des données, le message contient une ''enveloppe'' qui contient

plusieurs champs:

- Source: rang du process qui envoie.

- Destination: rang du process qui reçoit
- Etiquette: entier qui identifie le message de manière unique.
- Communicateur: communicateur au sein duquel se fait
l'échange.
Les fonctions d'envoi et de reception décrivent de manière non-ambiguë
et portable la nature des données transmises.
Benoit Semelin 200

MPI_SEND et MPI_RECV:
communication de base.
program EX_1 ! (prog de 2 process)
Nombre
implicit none Début
de variables
du tableau d'envoi Type MPI
include ''mpif.h'' à envoyer
des variables
(adresse)
integer :: rang,err,tag1,tag2 à envoyer
integer, dimension( MPI_STATUS_SIZE ) :: statut
Numéro
real, DIMENSION(10) :: x,x_remote du process
destinataire.
call MPI_INIT (err)
tag1=1
tag2=2
Etiquette
call MPI_COMM_RANK (MPI_COMM_WORLD ,rang,err) du message
x = rand()
if(rang = = 0) then
call MPI_SEND(x,10,MPI_REAL,1,tag1,MPI_COMM_WORLD,err)
Communicateur
call MPI_RECV(x_remote,10,MPI_REAL,1,tag2,MPI_COMM_WORLD,statut,err)
endif
if(rang = = 1) then
call MPI_RECV(x_remote,10,MPI_REAL,0,tag1,MPI_COMM_WORLD,statut,err)
call MPI_SEND(x,10,MPI_REAL,0,tag2,MPI_COMM_WOLD,err)
endif
Statut de
print*,'rang: ',rang, ' liste: ',x,x_remote l'opération de
call MPI_FINALIZE (err) Numéro réception.
du process
end program EX_0 source
Benoit Semelin 200

MPI_SEND et MPI_RECV:
détails d'utilisation.
- Pour des raisons de portabilité (grille hétérogène), il faut donner un ''type MPI''
pour les variables envoyées. Voici la correspondance pour Fortran:

INTEGER: MPI_INTEGER (MPI_INT en C)

REAL: MPI_REAL (MPI_FLOAT en C)
DOUBLE PRECISION: MPI_DOUBLE_PRECISION (MPI_DOUBLE en C)
LOGICAL: MPI_LOGICAL
CHARACTER(1): MPI_CHARACTER (MPI_CHAR en C)

- On peut aussi envoyer des types MPI_PACKED et des types définis par l'utilisateur
(structures).

- La taille du tableau (buffer) de reception doit être ≥ à la taille du message.

- Wildcards: MPI_RECV accepte MPI_ANY_SOURCE et MPI_ANY_TAG

(constantes mpi) comme valeur pour la source et l'étiquette du message.

- La variable ''statut'' contient (en Fortran):

- statut(MPI_SOURCE): source du message reçu.
- statut(MPI_TAG): l'étiquette du message reçu.
Benoit Semelin 200

Modes de communication:
définitions
Communication bloquante: l'appel à la fonction ne ''retourne'' que quand les ressources utilisées (p. e.
emplacement mémoire de la variable envoyée) peuvent être réutilisées et que la fonction ''complète''.
Communication non-bloquante: l'appel à la fonction retourne avant que les ressources aient été
libérées, et le programme continue. Il faut s'assurer qu'on ne modifie pas les ressources avant que la
communication soit effectivement complétée. La fonction ne complète qu'à ce moment là.

Mode d'envoi ''buffered'': le message est stoké dans une mémoire système locale avant d'être envoyé.
La fonction d'envoi (bloquante ou non) complète quand la copie est finie mais avant l'envoi.
Mode d'envoi ''synchrone'': l'envoi effectif du message ne commence que quand le process reçoit le
message « prêt » d'une commande de réception correspondante. La fonction d'envoi ne complète qu'à ce
moment là.
Mode d'envoi ''ready'': le process émetteur suppose que le recepteur est près à recevoir sans vérifier, et
envoie le message (meilleures performances). Si le récepteur n'a pas exécuté la commande de réception
correcpondante -> Erreur! La fonction complète quand l'evoi est fini.

Mode d'envoi ''standard'': Suivant la disponibilité en mémoire système, MPI choisit lui même entre
les modes ''buffered'' et ''synchrone''.

Il n'y a qu'un mode de reception.

MPI_SEND est une communication bloquante en mode standard.
MPI_RECV est une communication bloquante.
Benoit Semelin 200

Autres fonctions de
communication point à point
Il existe une panoplie de fonctions pour effectuer des communications
bloquante ou non dans le différent mode. Elle s'utilisent avec les même
arguements que MPI_SEND et MPI_RECV, plus un argument ''request''
pour les fonctions non-bloquantes.

Bloquant Non-bloquant
Standard MPI_SEND MPI_ISEND
Synchrone MPI_SSEND MPI_ISSEND
Ready MPI_RSEND MPI_IRSEND
Buffered MPI_BSEND MPI_IBSEND (nécessite MPI_BUFFER_ATTACH)

- Le mode buffered nécessite une copie mémoire de plus mais permet de continuer les calculs à
coup sûr.
- Le mode ready diminue la latence, mais est délicat à utiliser (synchronisation).
Benoit Semelin 200

Gérer la complétion
d'une opération non-bloquante.
Syntaxe d'un envoi non-bloquant:
MPI_ISEND(val,count,datatype,dest,etiquette,comm,requete,err)
Le code peut continuer mais garde la trace de l'envoi grâce à requete (poignée vers
un objet MPI). Il peut utiliser requete ultérieurement pour contrôler si la
communication est complète:
MPI_WAIT(requete,statut,err): attentant que la communication associée à requete soit
complétée. statut contient des infos sur la communication.
MPI_TEST(requete,flag,statut,err): flag=true si la communication associée à requete est
complétée. Sinon, flag=false, et requete est désallouée! Pas de vérification ultérieure possible.
MPI_REQUEST_GET_STATUS(request,flag,status,err): teste la complétion sans
désallouer la requête, même si flag=true.
MPI_REQUEST_FREE(request,err): désalloue la requête.

Il est possible de gérer des complétions mulitples grâce à MPI_WAITANY,

MPI_WAITALL, MPI_WAITSOME et les équivalent pour TEST.
Attention: MPI_REQUEST_FREE n'anulle pas la communication non bloquante.
Pour cela il faut utiliser MPI_CANCEL(resquest,err).
Benoit Semelin 200

Communications
non-prédictibles
Dans certains algorithmes, le besoin de communications peut dépendre
des données initiales (ex: modif dynamique des domaines):
- Les besoins de l'envoi sont déterminés localement par le calcul.
- Comment déterminer les réceptions à effectuer?
On peut tester périodiquement l'arrivée de messages sans les recevoir.
MPI_IPROBE(source,etiquette,comm,flag,status): non bloquant!
flag détermine si un message est arrivé ou non. source et etiquette
peuvent prendre les valeurs MPI_ANY_SOURCE et MPI_ANY_TAG,
si flag=true, status contient la source et l'étiquette.

Et si on ne connait pas la longueur du message?

- L'étiquette peut servir à coder le nombre et le type de donnée.
- Une précommunication à un format standardisé peut annoncer l'étiquette, le type et
le nombre de données d'un message à venir.
- Il est peut-être temps de passer à des communications « one-sided ».
Benoit Semelin 200

Synchronisation

Synchronisation globale:

MPI_BARRIER(comm,ierr): fonction collective

Permet de bloquer les process du communicateur comm jusqu`à ce que

le dernier soit arrivé à la barrière.

Synchronisation locale:
On peut utiliser une communication synchrone bloquante
(MPI_SSEND/RECV) pour synchroniser 2 process.
Si on a besoin de synchroniser de façon répétée un sous-groupe de
process, il faut sans-doute définir un nouveau sous-groupe-MPI et un
nouvel intra-communicateur.
Benoit Semelin 200

MPI:

communications globales
(fonctions collectives)
Benoit Semelin 200

Communication de
type ''broadcast''
Il s'agit de diffuser aux autres
processeurs une valeur connue sur
un seul:
MPI_BCAST(address,count,datatype,root,comm,err) CPU 1
X0
root désigne le rang du process
qui diffuse l'information. Les
autres arguments ont la même
signification que dans
MPI_SEND. CPU 0 CPU 2
X0 X0
C'est une fonction collective :
elle doit être appelée par tous les
process.

Il faut privilégier les

communications collectives pour
CPU 3
la simplicité de la programmation. X0
Benoit Semelin 200

Communication de
type ''gather''
Il s'agit de rassembler sur l'un des process
des données réparties sur l'emsemble des
process du communicateur :
MPI_GATHER(s_add , s_count , s_type , r_add , r_count , r_type , root , comm , err)

- root désigne le rang du process qui reçoit

les données. r_add désigne l'adresse (nom
CPU 1
de variable en fortran) où les données reçues X1
sont stockées.
- Les données sont stockées dans l'ordre de
rang des processeurs qui envoient.
X0
- r_count est le nombre de variables de type
CPU 0
X1
CPU 2
r_type reçue de chaque process. Donc, la
plupart du temps: X2 X2
s_count=r_count=constante ! X3

- Sur le process root, s_add peut être

remplacé par MPI_IN_PLACE. On suppose CPU 3
alors que les données à envoyer pour le
process root sont déjà à leur place dans le X3
tableau de reception.
Benoit Semelin 200

Autres communications de
type ''gather''
Pour recevoir une quantité de données différente de chaque process ou les disposer de
manière non-consécutive dans le tableau de réception, on utilise la version vecteur:

MPI_GATHERV(s_add , s_count , s_type , r_add , r_counts , disp , r_type , root , comm , err)

- s_count peut maitenant être différent sur chaque process.

- r_counts est un tableau de NB_PROCS entiers.
- r_counts(i) doit avoir la valeur de s_count sur le process i.
- disp est un tableau de NB_PROCS entiers.
- Les données reçues du process i sont stockée à l'adresse: r_add+disp(i)*sizeof(r_type)

Pour réaliser une opération sur les données reçues, on utilise une fonction de réduction:

MPI_REDUCE(s_add , r_add , count , datatype , op , root , comm , err)

- Les valeurs dans s_add sur les différents process sont combinées, élément à élément si count≠1,
et stockées dans r_add sur le process root. La combinaison est faite par l'opérateur op.
- En fortran op peut être: MPI_SUM, MPI_PROD, MPI_MAX, MPI_MIN, etc...
- Il est possible de créer ses propres opérateur grâce à MPI_OP_CREATE.
Benoit Semelin 200

Communications de
type ''scatter''
Il s'agit de répartir sur les process du
communicateur des données présentes sur
l'un des process:
MPI_SCATTER(s_add , s_count , s_type , r_add , r_count , r_type , root , comm , err)
- root désigne le rang du process qui
envoie les données. s_add désigne CPU 1
l'adresse (nom de variable en fortran) où
les données à envoyer sont stockées. X1
- Les données sont envoyées par paquet de
s_count aux process du communicateur
par ordre de rang.
X0
- En général r_count = s_count. CPU 0 CPU 2
X1
- Il existe une variante ''vecteur'':
X2 X2
MPI_SCATTERV
X3
- Il existe aussi:
MPI_REDUCE_SCATTER
Opére une reduction sur des tableaux, CPU 3
élément à élément, et stocke le résultat
pour les i-ème éléments sur le process de X3
rang i.
Benoit Semelin 200

Communications de
type ''allgather''
Même fonction que MPI_GATHER,
mais chaque process reçoit le résultat:
MPI_ALLGATHER(s_add , s_count , s_type , r_add , r_count , r_type , comm , err)

- Pas de root !
CPU 1
- Autres arguments identiques à
X0 X1 X2X3
MPI_GATHER.
- Equivalent à NB_PROCS appels à
MPI_GATHER avec à chaque fois un
process différent comme root. X0 X0

CPU 2
CPU 0
- Il existe une version vecteur: X1 X1
MPI_ALLGATHERV X2 X2
- Et une version avec réduction: X3 X3

MPI_ALLREDUCE
X0 X1 X2 X3

CPU 3
Benoit Semelin 200

Communications de
type ''all-to-all''
Même fonction que
MPI_ALLGATHER, mais chaque
process reçoit des données
différentes:
CPU 1
MPI_ALLTOALL(s_add , s_count , s_type ,
r_add , r_count , r_type , comm , err) W1 X1 Y1 Z1
X0 X1 X2 X3
- A rg u m e n t s i d e n t i q u e s à
MPI_ALLGATHER.

- Mais ici s_add et r_add

désignent des tableaux de même CPU 0 W0W0 Y0W2
taille.

CPU 2
X0 W1 Y1 X2
- Il existe deux versions vecteurs, Y0 W
2
Y2 Y2
pour pouvoir faire varier le
Z0 W3 Y3 Z
nombre, le type et la position des 2
données envoyées et reçues:
MPI_ALLTOALLV Z0 Z1 Z2 Z3
MPI_ALLTOALLW
W3 X3 Y3 Z3
Avec MPI_ALLTOALLW on peut
tout faire! Ou presque... CPU 3
Benoit Semelin 200

Types dérivés:

variables MPI
définies par l'utilisateur
Benoit Semelin 200

Déclarer un type dérivé

Les fonctions de communication MPI utilisent comme arguments des types MPI
standard comme MPI_INTEGER ou MPI_REAL. Il est possible de définir des types
plus complexes. Cela peut servir, par exemple:
- A envoyer une section de tableau.
- A envoyer une variable de type structure définie dans le programme.
Exemple de déclaration d'un type dérivé simple:

INTEGER :: err,MPI_vector
call MPI_TYPE_CONTIGUOUS(3,MPI_REAL,MPI_vector,err)
call MPI_TYPE_COMMIT(MPI_vector,err)
... communications ...
call MPI_TYPE_FREE(MPI_vector,err)
On définit ici un type MPI_vector, constitué de 3 MPI_REAL stockés consécutivement
en mémoire. MPI_vector peut servir à définir d'autres types dérivés.
Il faut ''compiler'' le type avec MPI_TYPE_COMMIT avant de pouvoir l'utiliser dans
une communication.
Benoit Semelin 200

Transmettre une section de

tableau
Les sections de tableaux sont un outils puissant de fortran. On peut définir des types dérivés pour
envoyer des sections de tableau. Exemple:

Les éléments foncés constituent la section x(2:12:3) du tableau x(1:12). Définissons un type MPI
correspondant.
INTEGER :: err,MPI_vector_section
REAL, DIMENSION(12) :: x,y

call MPI_TYPE_VECTOR(4,1,3,MPI_REAL,MPI_vector_section,err)
call MPI_TYPE_COMMIT(MPI_vector_section,err)

call MPI_ALLREDUCE(x(2),y(2),1,MPI_vector_section,MPI_SUM,MPI_COMM_WORLD,err)

MPI_TYPE_VECTOR construit le type MPI_vector_section, constitué de 4 blocs de 1 MPI_REAL,

avec un pas entre les débuts de blocs de 3 MPI_REAL. 4, 1 et 3 peuvent varier à volonté.
Attention: les tableaux multidimensionnels sont stockés sous forme de vecteur 1D. Il faut savoir quel est
l'indice qui varie le plus vite. On peut alors définir des types sections recursivement sur les dimensions.

On peut définir des sections avec pas variables entre blocs avec MPI_TYPE_CREATE_INDEXED_BLOCK,
et avec pas et longueusr de blocs variables avec MPI_TYPE_INDEXED.
Benoit Semelin 200

Type dérivé MPI

correspondant à une struture
integer :: err,MPI_integer_length,MPI_real_length,MPI_logical_length
integer :: MPI_vector_length,MPI_vector,MPI_particle
integer, dimension(10) :: array_of_block_length,array_of_types,array_of_displacement

call MPI_TYPE_EXTENT(MPI_INTEGER,MPI_integer_length,err)
call MPI_TYPE_EXTENT(MPI_REAL,MPI_real_length,err)
call MPI_TYPE_EXTENT(MPI_LOGICAL,MPI_logical_length,err)
} Pour la portabilité
(taille mémoire en bits)

}
call MPI_TYPE_CONTIGUOUS(3,MPI_REAL,MPI_vector,err)
call MPI_TYPE_COMMIT(MPI_vector,err) Définition d'un type vecteur
call MPI_TYPE_EXTENT(MPI_vector,MPI_vector_length,err)

}
array_of_block_length(1:3) = (/1,2,1/)
array_of_types(1:3) = (/MPI_INTEGER,MPI_vector,MPI_REAL/)
array_of_displacement(1) = 0 Description du
array_of_displacement(2) = MPI_integer_length type structure
array_of_displacement(3) = array_of_displacement(2) + 2*MPI_vector_length (déplacements en bits)

call MPI_TYPE_STRUCT(3,array_of_block_length(1:3),array_of_displacement(1:3) &

& ,array_of_types(1:3),MPI_particle,err)
call MPI_TYPE_COMMIT(MPI_particle,err)
} Déclaration
du type

Pas de moyen de définir des types objets dans MPI de base!

Benoit Semelin 200

MPI:

communications ''one-sided''
Benoît Semelin 200

Introduction

Que faire quand le process receveur ne sait pas qu'il doit recevoir?

Lancer régulièrement de communications globales.

Performances?
Utiliser MPI_IPROBE et des communications non bloquantes...
synchronisation délicate.
Utiliser des communications ''one-sided''.

Les communications one-sided permettent d'utiliser le caractère SMP

d'une architecture pour améliorer les performances, mais ne constituent
pas un modèle complet de programmation SMP.
Benoit Semelin 200

Déclarer une « fenêtre »

de communication one-sided.
Il faut définir une zone mémoire accessible par les process distants lors
de communications one-sided:
MPI_WIN_CREATE(add , size , disp_unit , info , comm , win , err)
Ø add: addresse de début de la fenêtre (ex, nom de tableau).
Ø size: taille de la fenêtre en bit.
Ø disp_unit: unité de déplacement pour les accès ultérieurs.
Ø info: poignée vers un tableau où les 3 premiers seront stockés.
Ø comm: communicateur
Ø win: poignée (entier) attachée à la fenètre créée.
Ø err: code d'erreur
MPI_WIN_CREATE est une fonction collective. La fenêtre n'a pas
forcement la même taille et position sur tous les process.
MPI_WIN_FREE(win,err) pour désallouer la fenêtre.
Benoit Semelin 200

MPI_PUT / MPI_GET:
écriture / lecture
Il y a 2 fonctions principales pour réaliser une communication one-sided:
MPI_PUT(loc_add , loc_count , loc_type , remote_rank , remote_disp &
& , remote_count , remote_type , win , err)
MPI_GET(loc_add , loc_count , loc_type , remote_rank , remote_disp &
& , remote_count , remote_type , win , err)
Ø loc_add: adresse de début des données sur le process qui appelle la fonction
Ø remote_disp: Définit l'endroit où lire/écrire les données par le déplacement en
disp_unit par rapport au début de la fenêtre sur le process cible.
Ø loc_count,remote_count: Nombres de variables à lire/écrire.
Ø loc_type,remote_type: Type des variables à lire/écrire.
Ø remote_rank: Numéro du process cible.
Ø win: fenêtre dans laquelle se fait la communication.
Il s'agit de communications non bloquantes ! Il faut vérifier la complétion par des
appel de synchronisation.
Il existe aussi MPI_ACCUMULATE.
Benoit Semelin 200

Complétion des communications

one-sided.
CALL MPI_WIN_CREATE(x , 1000, MPI_real_length, info , MY_COMM_WORLD, x_win , err)
CALL MPI_WIN_FENCE(MPI_MODE_NOPRECEDE,x_win,err)

target=mod(rank+1,nb_proc)
index_to_put=local_computation()
CALL MPI_PUT(x(index_to_put),1,MPI_REAL,target,index_to_put,1,MPI_REAL,x_win,err)

CALL MPI_WIN_FENCE(MPI_MODE_NOSUCCEED,x_win,err)
CALL MPI_WIN_FREE(x_win,err)

MPI_WIN_FENCE est une fonction collective. Elle agit comme une barrière, aucun process
du groupe ne continue tant que toutes les communications de la fenêtre ne sont pas complétées.
Le premier argument sert à l'optimisation, il peut toujours valoir 0. Il peut prendre les valeurs:
- MPI_MODE_NOSTORE: pas d'écriture en local dans la fenêtre depuis le dernier FENCE.
- MPI_MODE_NOPUT: Pas de PUT vers la fenêtre locale d'ici le prochain FENCE.
- MPI_MODE_NOPRECEDE: Pas de communications antérieures à compléter.
- MPI_MODE_NOSUCCED: Pas de communications ultérieures.
On peut combiner ces valeurs (voir manuel de réference). Les optimisations correspondantes ne
sont parfois pas implémentées.
Il est possible de synchroniser les process 2 à 2 avec MPI_WIN_START,
MPI_WIN_COMPLETE, MPI_WIN_POST et MPI_WIN_WAIT.
Benoit Semelin 200

Ce que je n'ai pas traité:

- Topologie de process
- Entrés-sorties parallèles
- Inter-communications
- Création dynamique de process
- Packing de données
 GPU

Global Memory
PC I/O: 4 GB.s-1
Shared mem

Slow
Fast
Shared mem Shared mem
CU Local mem, reg

CU Local mem, reg

Multi CU Local mem, reg Multi
processor processor
Multi processor

Exemple: Nvidia tesla C1060

- Global memory: 4 GB - 30 mulprocessors
- Shared memory: 16 kB/multiproc ! - 8 CU per multiprocessor
Ø CUDA: extension de language
§ C : nvcc (Nvidia)
§ Fortran: pgf90 (depuis Nov 2009)

Ø OpenCL: librairies.
Tout neuf! 1 implementation?

Ø Directives de compilation (style OpenMP)

Compilateurs PGI (pgf90, pgcc).
ü Ecrire un code scalaire pour CPU

ü Identifier les subroutines coûteuses

ü Vérifier le caractère presque Single Instruction Multiple Data

IF-THEN-ELSE est acceptable mais « sérialisé »

ü Copier les données dans la mémoire globale du GPU

ü Envoyer la subroutine pour exécution sur le GPU

ü Récupérer les résultats depuis la mémoire du GPU.

Sur le GPU: N treads exécutent chacune toutes les instructions du kernel.

- Les threads sont organisées en blocs.

Threads d’un bloc:
- accès à une mémoire partagée commune (16 ko)
- organisées en grille 1D, 2D ou 3D.

- thread identifiée par: threadIdx.x, threadIdx.y, threadIdx.z (en C)

threadIdx%x, threadIdx%y, threadIdx%z (en Fortran)
- Nb max de thread par bloc: 512 (dim max en z 64)

- On définit une grille de block ( < 655353 ! mais si 30 multiprocesseurs, 30 blocs actifs)
Bloc identifié par: blockIdx.x, blockIdx.y, blockIdx.z (en C)
blockIdx%x, blockIdx%y, blockIdx%z (en Fortran)

threadIdx et blockIdx servent à paramétrer l’exécution.

 #include <cuda_runtime.h> MODULE GPU_KERNELS
#include <cude_runtime_api.h> USE CUDAFOR
CONTAINS
ATTRIBUTE(global) SUBROUTINE ma_subroutine_kernel(args);
__global__ void ma_function_kernel(args);
integer :: i
i=griddim%x+blockdim%y-threadidx%x/blockidx%z
void main( ) { END SUBROUTINE ma_subroutine_kernel
END MODULE GPU_KERNEL
dim3 dimBlock(512); !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
dim3 dimGrid(4,4,4); PROGRAM TOTO
USE CUDAFOR
USE GPU_KERNEL
ma_fonction_kernel<<<dimGrid,dimBlock>>>(args);
} type(dim3) :: dimblock, dimGrid

__global__ void ma_function_kernel(args); dimGrid=dim3(4,4,4)

{ Dimblock=dim3(512,1,1)
int i; call ma_subroutine_kernel<<<dimGrid,dimBlock>>>(args)
i=gridDim.x+blockDim.y-threadIdx.x/blockIdx.z;
} END PROGRAM

Depuis un kernel, gridDim, blockDim, threadIdx et blockIdx sont prédéfinies.

Types possible de qualificateurs de functions/subroutines:
host (defaut, execution sur CPU)
global(GPU, appelable depuis CPU seulement)
device (GPU, appelable depuis GPU seulement).
#include <cuda_runtime.h> MODULE GPU_KERNELS
#include <cude_runtime_api.h> USE CUDAFOR
#define NCELL 8192 CONTAINS
ATTRIBUTE(global) SUBROUTINE ma_subroutine_kernel(vel,ncell);
__global__ void ma_function_kernel(args); real :: vel (0,ncell-1)
integer, value :: ncell
void main( ) { vel(dimGrid%x*BlockIdx%x+threadIdx%x)= …..
END SUBROUTINE ma_subroutine_kernel
float v[NCELL]; END MODULE GPU_KERNEL
float* v_dev; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
PROGRAM TOTO
dim3 dimBlock(512,1,1); USE CUDAFOR
dim3 dimGrid(32,1,1); USE GPU_KERNEL
integer, parameter :: ncell=8192
cudaMalloc( (void**)&v_dev , ncell*sizeof(float)); type(dim3) :: dimblock, dimGrid
cudaMemcpy(v_dev, &v, ncell*sizeof(float), cudaMemcpyHostToDevice);
ma_fonction_kernel<<<dimGrid,dimBlock>>>(v_dev); real(KIND=4), dimension(0:ncell-1) :: v
cudaMemcpy(&v,v_dev,ncell*sizeof(float),cudaMemcpyDeviceToHost); real(KIND=4), device, dimension(0:ncell-1) :: v_dev
} v=…
v_dev=v
__global__ void ma_function_kernel(float* vel_in_kernel);
dimGrid=dim3(32,1,1)
{
Dimblock=dim3(512,1,1)
int i;
vel_in_kernel[dimGrid.x*BlockIdx.x+threadIdx.x]= …..;
call ma_subroutine_kernel<<<dimGrid,dimBlock>>>(v_dev,ncell)
}
v=v_dev
END PROGRAM
v_dev est sur la mémoire globale du GPU.

Variables déclarés dans un kernel: mémoire locale du multiprocesseur (valeur ≠ pour chaque thread)
 Accès à la mémoire globale: ~ qq centaines de cycles d’horloge
Accès à la mémoire partagé: ~1 cycle d’horloge

⇒Si plusieurs accés à une variable, la copier en mémoire partagée.

#include <cuda_runtime.h> MODULE GPU_KERNELS
#include <cude_runtime_api.h> USE CUDAFOR
#define NCELL 8192 CONTAINS
ATTRIBUTE(global) SUBROUTINE ma_subroutine_kernel(vel,ncell);
__global__ void ma_function_kernel(args); real :: vel (0,ncell-1)
integer, value :: ncell
void main( ) { real, shared :: vel_sh(0:511)
….. vel_sh(threadIdx%x)=vel(dimGrid%x*BlockIdx%x+threadId%x)
} …..
END SUBROUTINE ma_subroutine_kernel
__global__ void ma_function_kernel(float* vel_in_kernel); END MODULE GPU_KERNEL
{ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
__shared__ float vel_sh[512]; PROGRAM TOTO
vel_sh[threadIdx.x]=vel_in_kernel[dimGrid.x*BlockIdx.x+threadIdx.x]; …..
}
END PROGRAM

Tenter d’effectuer des accès mémoire coalescents !

Pour synchroniser les threads d’un bloc (dans un kernel):
En C: __syncthreads()
En fortran: call syncthreads()

Pour synchroniser les blocs à la sortie d’un kernel :

En C: cudaThreadSynchronize()
En fortran: err=cudaThreadSynchronize()

(avant de lire dans la mémoire globale du GPU?)

Equation de Burger: turbulence 1D ∂v ∂v 1 ∂ 2v
+v =
∂t ∂x Re ∂x 2

1 solution particulière: Conditions numériques:

Ø 8192 cellules
Ø 2 106 pas de temps

Resultat du Benchmark:
ü gcc –O3: 461 s
ü icc –O3: 227 s
ü pgf90 –fast: 309 s
ü ifort –O3: 38.8 s
ü nvcc, pgf90 -Mcuda : 3.67 s
- Les performances dépendent de:
§ Haute intensité arithmétique

§ Accés mémoire coalescent

§ Maximisation de la bande passante mémoire (mem partagée)

- Plus facile avec les codes sur grille.

- Il existe une librairie cuda standard:

FTT, tri, histogramme, générateur aléatoire, etc…

- Ray-tracing for gravitational lensing (speed up 100)

- Cosmological interacting scalar field (speed up 10 to 100)

- Adaptive grid hydro (speed up ~ 10)

- Digital correlator for radio-interferometer (MWA) (speed up 60)

- MHD (speed up 43)

- Direct N-body and tree code (speed up > 10)

- Radiative transfer (moment method, D. Aubert) (speed up 80)