Document Actions
EDP parallélisées
Exemples de parallélisation d'EDP par OpenMP et MPI
Introduction
Cette page contient quelques exemples de parallélisation d'Équations aux Dérivées Partielles (EDP) simples.Le but est de montrer comment on utilise MPI et OpenMP pour partager les calculs entre différents processus sur différents processeurs.
Les exemples traités sont volontairement simples qu'il s'agisse des équations, des algorithmes ou des schémas numériques employés.
Une page particulière est consacrée à chaque exemple :
- définition du problème traité,
- cas test et validation,
- mesures de performances du code,
- mesures de l'efficacité de la parallélisation,
- accès aux codes sources et autres fichiers de résultats.
Choix des algorithmes
On est souvent amené à résoudre des systèmes linéaires (possiblement creux) de dimension assez large au sein de problèmes complexes. Dans le cas de programmes parallèles, le choix de méthodes de résolution de systèmes linéaires doit aussi prendre en compte le potentiel de parallélisation de ces algorithmes. Étant donnés les problèmes de dépendance des factorisations et des résolutions de systèmes triangulaires, les méthodes directes n'ont pratiquement pas de parties parallélisables. Les méthodes itératives ont des possibilités naturelles de parallélisme plus grandes. Elles reposent sur le produit matrice vecteur qui, en général, se prête bien au parallélisme.C'est pour cela que les méthodes choisies ici sont deux méthodes itératives, dont on ne présente que les algorithmes.
Méthode de Jacobi
Il s'agit d'une méthode très simple à programmer pour résoudre un système linéaire de la forme A.x = b.Le stockage en mémoire est assez limité et lorsqu'elle converge il faut que celle-ci soit suffisamment rapide pour que le nombre d'itérations ne soit pas trop grand.
Pour A = (a_ij), b = (b_i) et x^{(0)} = (x_i), on construit une suite de vecteurs x^{(k)} à l'aide de la relation
Le calcul de l'itéré (k+1) en fonction de l'itéré (k) est proche du produit matriciel de l'itéré (k) par A en terme de dépendances des données et du nombre d'opérations et pourra donc être parallélisé.
Le critère d'arrêt de la méthode est sur
il est nécessaire de s'assurer que la méthode converge avant de l'appliquer à un système linéaire. Par exemple, si la matrice A est à diagonale strictement dominante alors la méthode de Jacobi est convergente.
Méthode du Gradient conjugué
La méthode du Gradient Conjugué est une méthode de descente qui permet d'obtenir une solution approchée d'un système linéaire en minimisant une fonctionnelle. La matrice A du système à résoudre est supposée symétrique et définie positive.En notant (a,b) le produit scalaire euclidien et ||x|| la norme associée, l'algorithme est le suivant :
Soit le vecteur x_0 donné et p_0 = r_0 = b - A.x_0.
Pour k = 0, 1, ...
Le critère d'arrêt porte sur
C'est un peu plus complexe que la méthode de Jacobi mais beaucoup plus efficace dans le cas d'une matrice symétrique définie positive.
Une itération du Gradient Conjugué est constituée d'un produit matrice-vecteur, de deux produits scalaires et de trois saxpy (combinaison linéaire de vecteurs). Toutes ces opérations sont parallélisables.
Remarque :
il est nécessaire de s'assurer que la matrice A est symétrique, définie positive pour appliquer la méthode du Gradient Conjugué. Cette propriété est vérifiée dans nos exemples.
Choix des schémas numériques
Les schémas numériques employés pour la discrétisation spatiale et pour la discrétisation temporelle sont simples à mettre en oeuvre et permettent de se concentrer sur les techniques de parallélisation. Ils peuvent, bien entendu, être remplacés par des schémas plus performants.Discrétisation en espace
Le schéma spatial employé est de type différences finies d'ordre deux à trois points. Le cas qui nous intéresse est la discrétisation du laplacien en 3 dimensions d'espace :
Sur un maillage uniforme noté x_{i,j,k} avec des points espacés du pas de discrétisation hx, hy, hz respectivement selon les directions spatiales, cela s'écrit :
Ce schéma est consistant à l'ordre 2.
Discrétisation en temps
Nous allons considérer successivement plusieurs problèmes : les deux premiers sont stationnaires et le troisième est instationnaire.Seul le troisième est discrétisé en temps et à l'aide du schéma d'ordre deux semi-implicite de Crank-Nicholson. La discrétisation en temps de l'équation
s'écrit
Dans le cas qui nous intéresse, ce schéma est inconditionnellement stable.
Globalement, nous avons une discrétisation temps - espace consistante à l'ordre 2 et stable, donc convergente.
Choix des parallélisations
Deux approches sont considérées dans les exemples traités :Décomposition de domaines
Cette technique de parallélisation consiste à partager le domaine de calcul entre plusieurs morceaux, appelés sous-domaines, dont la somme directe forme le domaine initial. Il s'agit d'une parallélisation explicite puisque l'ensemble du travail reste à la charge du programmeur. Cette technique s'applique aussi bien sur des maillages structurés que non structurés.Que la discrétisation spatiale soit faite par Volumes ou Éléments Finis (i.e. cellules de calcul) ou Différences Finies (noeuds de calcul), il est très souvent nécessaire de connaitre des valeurs autour de la maille ou du noeud de calcul courant pour déterminer les valeurs dans la maille ou sur le noeud courant (autour du domaine initial, ces valeurs sont données par les conditions aux limites).
Lors de la décomposition du domaine, chaque sous-domaine reçoit donc en plus de ses cellules (ou noeuds) une ou plusieurs couches de cellules (ou noeuds) fantomes provenant de ses sous-domaines voisins et permettant le calcul des valeurs pour ses propres cellules situées sur son bord. Ces cellules sont mises à jour à chaque pas de temps avec les valeurs calculées par le processus qui traite le sous-domaine propriétaire de ces cellules.
Cette technique est souvent appliquée avec MPI (Message Passing Interface) pour gérer les échanges de valeurs aux frontières entre sous-domaines adjacents avec une topologie. Cela peut aussi se faire avec OpenMP mais c'est plus technique à concevoir.
Prenons comme exemple, un domaine 2D rectangulaire avec un maillage uniforme.
On découpe le maillage en sous-domaines rectangulaires de manière classique et chaque processus MPI reçoit un tel bloc. On superpose une grille de processus sur la grille des sous-domaines grace à une topologie construite par MPI (MPI_DIMS_CREATE, MPI_CART_CREATE). La bibliothèque MPI possède des fonctions permettant de connaitre les processus voisins (MPI_CART_SHIFT) et de se repérer dans cette décomposition (MPI_CART_COORDS et MPI_CART_RANK).
Ainsi, chaque sous-domaine possède un ou plusieurs voisins et donc des cellules fantomes provenant de ces voisins le long de leurs frontières communes. Pour chacune de ces interfaces a lieu un échange de données et la mise à jour est faite au moyen d'envois/réceptions standards (MPI_SEND, MPI_RECV).
Dans le cas représenté, chaque sous-domaine est un rectangle et les interfaces sont situées le long de ses cotés. Pour une interface donnée, chaque processus situé de part et d'autre envoie et reçoit des données; on utilise alors la fonction MPI_SENDRECV. Dans le cas de données stockées de manière discontigue, il est parfois judicieux de construire un type dérivé pour réaliser les communications.
Dans nos exemples, le schéma spatial fait uniquement intervenir les points directement adjacents au noeud de calcul courant. Ainsi pour chaque sous-domaine il n'y a qu'une couche de noeuds fantomes.
Les conventions de notations sont alors les suivantes, par exemple pour la direction d'espace x (de même pour les directions y et z) :
- indices globaux de calcul du domaine : de 1 à ntx;
- indices locaux de calcul pour chaque sous-domaine : de sx à ex;
- indices locaux de stockage pour chaque sous-domaine : de (sx-1) à (ex+1) (NB : longueur = ex-sx+3).
- Pour le plan xOy, il s'agit de (ey-sy+1) vecteurs régulièrement espacés (longueur d'un vecteur : ex-sx+1, pas entre deux débuts consécutifs : ex-sx+3), donc on crée le type ip_xy avec la fonction MPI_TYPE_VECTOR.
- Pour le plan xOz, il s'agit de (ez-sz+1) vecteurs régulièrement espacés (longueur d'un vecteur : ex-sx+1, pas entre deux débuts consécutifs : (ex-sx+3)*(ey-sy+3) ), donc on crée le type ip_xz avec la fonction MPI_TYPE_VECTOR.
- Pour le plan yOz, il s'agit d'éléments isolés irrégulièrement espacés, en fait deux pas interviennent. Un pas constant dans la direction y et un pas constant dans la direction z. On crée deux types dérivés, le second basé sur le premier en mesurant les distances entre les éléments en octets à l'aide de la fonction MPI_TYPE_HVECTOR.
|
|
| CALL MPI_TYPE_SIZE (ITYPE_REAL, taille_real, ierr ) |
|
|||||||||||||
|
Le calcul d'un saxpy ne fait intervenir que les valeurs associées aux noeuds de calcul, donc pas celles des noeuds fantomes et par conséquent il ne nécessite pas de mise à jour des valeurs par des communications.
Le calcul d'un produit scalaire est distribué entre les processus, chacun d'eux réalisant le produit sur la portion de vecteur correspondant à son sous-domaine. Une opération de réduction impliquant l'ensemble des processus permet d'obtenir la valeur globale du produit scalaire, résultat ensuite transmis à tous les processus au travers d'une diffusion collective (MPI_ALLREDUCE).
Pour le produit matrice vecteur (ou assimilé) les calculs font intervenir le voisinage des noeuds de calcul. A l'intérieur des sous-domaines, le voisinage est connu puisque dans le sous-domaine lui-même. Par contre pour les noeuds situés au bord, le produit fait intervenir les valeurs des noeuds fantomes et cela nécessite une phase de communications entre tous les sous-domaines pour mettre à jour les valeurs de ces noeuds.
Remarque :
Cette technique de parallélisation peut aussi être appliquée avec OpenMP.
Pour cela il suffit de créer à la main la décomposition en sous-domaines, de gérer les "échanges" de données aux interfaces de sous-domaines. Pour la première partie on peut s'inspirer de ce qui est fait avec les fonctions MPI de gestion de topologies. Pour la seconde partie, on crée des tableaux buffers qui contiennent les valeurs à transmettre aux voisins et on travaille avec des variables de type POINTER qui contiennent les adresses mémoire de ces buffers. De là, il suffit d'aller lire les données stockées dans les buffers des voisins correspondants.
Dans ce cas, la parallélisation est explicite et on ne met plus de directives de la forme !$OMP DO autour des boucles puisque chaque thread OpenMP gère son propre sous-domaine. Seules les directives de synchronisation, les barrières, les mises à jour atomiques persistent.
Partage du travail
Dans le cas de boucles de calcul, cette stratégie repose sur la distribution des itérations entre les processus. Il s'agit d'une technique de parallélisation implicite puisque les synchronisations, distribution effective des itérations sont réalisées par la bibliothèque parallèle. Un exemple est OpenMP, ensemble de directives de parallélisation pour architecture à mémoire partagée. Tout le travail de la parallélisation se trouve indiqué au compilateur au travers de directives insérées dans le code source. Elles sont préfixées par une sentinelle, de la forme !$OMP ou C$OMP.L'aspect important est la visibilité des variables :
- une variable est publique (SHARED) si elle est commune à tous les processus ;
- une variable est privée (PRIVATE) si chaque processus dispose de sa copie privée.
Il faut distinguer les variables statiques (emplacement en mémoire défini dès sa déclaration) des variables automatiques (emplacement attribué lors du lancement de l'unité de programme contenant leur déclaration jusqu'à la fin de son exécution); les variables globales (définies dans le programme principal) des variables locales. Une variable locale initialisée à sa déclaration (valeur explicite, par DATA ou SAVE) est statique sinon elle est automatique. Il en est de même pour les tableaux : les tableaux automatiques sont des tableaux locaux dont les dimensions dépendent des arguments reçus alors que les dimensions des tableaux dynamiques peuvent varier d'une exécution à l'autre et ils sont alloués à la demande explicite du programmeur.
De manière générale, il faut se souvenir que l'attribut implicite des variables d'un sous-programme appelé dans une région parallèle est :
SHARED pour les variables locales statiques, déclarées en common ou en module, passées en argument;
PRIVATE pour les variables locales et les tableaux automatiques.
Dans le cas d'un domaine 3D, le domaine n'est pas partagé mais les boucles de calcul, après analyse des dépendances des calculs, sont englobées dans des régions parallèles OpenMP et parallélisées par des directives de la forme !$OMP DO.
Un code parallélisé par OpenMP est un code séquentiel qui contient des régions parallèles. Chaque région voit la création d'une équipe de threads qui se partage le travail à effectuer au sein de cette région. La création d'une telle équipe n'est pas gratuite et il est souhaitable de réduire au minimum le nombre de régions parallèles. Dans notre cas, on va n'en créer qu'une seule et la placer le plus haut possible dans la hiérarchie des unités du programme.
La boucle de l'algorithme itératif ne peut pas être parallélisée puisqu'il y a une dépendance évidente entre deux itérations consécultives, elle sera donc dupliquée : chaque thread l'exécutant.
| !$OMP PARALLEL DEFAULT( NONE ) | |||||
| !$OMP SHARED(u,b,u_exact,f1,f2,p,q,r,rnrm2,rnrmoo) | |||||
| !$OMP SHARED(tps,dt,rnorm_k,itgc,nbitgc,rerr2,rerroo,prec) | |||||
| !$OMP SHARED(nbproc,nbiter,machep,ifreq,it_max,rnudt2,rlambda) | |||||
| !$OMP PRIVATE(itg) | |||||
| DO itg = 1, nbiter | |||||
| !$OMP MASTER | |||||
| nbitgc = nbitgc + itgc | |||||
| itgc = 0 | |||||
| rnorm_k = zero | |||||
| |||||
| CALL scdmb( b, f1, f2, u, tps, rnudt2 ) | |||||
| !$OMP MASTER | |||||
| tps = tps + dt | |||||
| |||||
| CALL gradconj( u, b, p, q, r, it_max, prec, itgc, rnorm_k, rnudt2 ) | |||||
| END DO | |||||
| !$OMP END PARALLEL |
La boucle du gradient conjugué est elle-aussi dupliquée : il y a dépendance évidente entre deux itérations consécultives :
| CALL pmv( r, u, - rnudt2 ) | ||||||||||||||||
| CALL saxpy( r, -one, b, one ) | ||||||||||||||||
| ||||||||||||||||
| convergence = ( SQRT(rnorm_k) < prec ) | ||||||||||||||||
| DO WHILE ( (.NOT. convergence) .AND. (it < it_max) ) | ||||||||||||||||
| it = it + 1 | ||||||||||||||||
| CALL pmv( q, p, - rnudt2 ) | ||||||||||||||||
| CALL prodscal ( alpha_k, q, p ) | ||||||||||||||||
|
||||||||||||||||
| CALL saxpy( u, one, p, alpha_k ) | ||||||||||||||||
| CALL saxpy( r, one, q, - alpha_k ) | ||||||||||||||||
| CALL prodscal( rnorm_k, r, r ) | ||||||||||||||||
|
||||||||||||||||
| CALL saxpy( p, beta_k, r, one ) | ||||||||||||||||
| convergence = ( SQRT(rnorm_k) < prec ) <-- convergence : variable locale, donc privée | ||||||||||||||||
| END DO |
Par contre, les boucles de calcul telles que les produits matrice vecteur, les produits scalaires (en premier ci-dessous), les saxpy (en second ci-dessous) sont parallélisées puisqu'il n'y a pas de dépendance.
produit scalaire :
| !$OMP DO REDUCTION( +:rnorm ) | |||
| DO k = sz, ez | |||
| DO j = sy, ey | |||
| DO i = sx, ex | |||
|
|||
| END DO | |||
| END DO | |||
| !$OMP END DO |
saxpy (combinaison linéaire de vecteurs) :
| !$OMP DO |
| DO k = sz, ez |
| DO j = sy, ey |
| DO i = sx, ex |
| a(i,j,k) = scala * a(i,j,k) + scalb * b(i,j,k) |
| END DO |
| END DO |
| END DO |
| !$OMP END DO |
Ainsi on retrouve des directives de la forme !$OMP DO dans des sous-programmes, dans des fichiers indépendants sans qu'elles soient comprises entre un !$OMP PARALLEL / !$OMP END PARALLEL, le tout dans un même sous-programme. On appelle cela des directives orphelines (orphan directives).
On distingue le produit scalaire des deux autres opérations puisque les données en entrée sont des vecteurs et la donnée en sortie est un scalaire, il y a donc une réduction. Avec OpenMP cela se traduit par l'ajout d'une clause de réduction au travers de la directive !$OMP REDUCTION (op:var) où op désigne l'opération effectuée (supposée associative et commutative) et var désigne la variable scalaire impliquée.
De même que l'on déclare les variables au début du programme (ou des sous-programmes), on indique la visibilité des variables lors de l'entrée dans la région parallèle OpenMP. Pour cela on utilise la clause DEFAULT(NONE) qui oblige le compilateur à s'assurer que toutes les variables ont été renseignées. Un mauvais choix peut entrainer des résultats erronés et différents d'une exécution à l'autre en parallèle alors que le code séquentiel est correct.
Les tableaux (solution calculée, second membre, vecteurs de travail) sont accessibles pour toutes les threads qu'ils soient lus ou écrits puisque les boucles qui les parcourent sont partagées. Par contre les indices de boucle, sont privés pour chaque thread, de même que certaines variables scalaires car c'est moins cher de faire faire un petit calcul par toutes les threads plutot que de mettre une synchronisation suite au calcul effectuer par une seule d'entre-elles.
De manière générale, les variables utilisées en lecture sont partagées (ou publiques, clause SHARED) alors que les variables locales écrites sont privées (clause PRIVATE).Il arrive que certaines variables scalaires soient modifiées et leur nouvelle valeur doit ensuite être accessible à l'ensemble des threads. On utilise alors des directives qui limitent leur accès en écriture à une seule thread au moyen des directives !$OMP MASTER / !$OMP END MASTER et !$OMP BARRIER qui synchronise les threads.
Cela intervient notamment dans l'algorithme du Gradient Conjugué pour les scalaires alpha_k et beta_k : une seule thread les calcule et une synchronisation les rend accessible lorsque la nouvelle valeur est disponible.
Validation
Une fois les programmes écrits, on procède par étapes de la manière suivante :- déboguage du code : options de déboguage
- optimisation du code : options d'optimisation
- calcul d'erreur avec solution analytique (voir exemples).
Comme la méthode de Jacobi converge très lentement, le critère d'arrêt fonctionne aussi sur un nombre maximal d'itérations.
Mesure de performances du code
Pour chaque algorithme, on évalue l'ordre de grandeur du nombre d'opérations effectuées en fonction :- des dimensions du problème : ntx, nty et ntz;
- du coût en nombre d'opérations d'un produit matrice vecteur, d'un produit scalaire, d'un saxpy;
- du nombre d'itérations de la méthode itérative;
- du nombre d'itérations en temps, le cas échéant.
Ces valeurs sont à comparer aux performances crête théorique des processeurs Power 4, 1300 MHz, 5200 MFlops des noeuds IBM.
Mesure d'efficacité de la parallélisation
Une technique simple consiste à considérer comme référence le temps écoulé T_s nécessaire pour la résolution du problème par un code séquentiel optimal (algorithmes et schémas numériques efficace, options d'optimisation). De là, on mesure le temps écoulé T(p), pour la résolution par le code parallèle avec p=2, 4, 8 processus.L'accélération (ou speedup) est définie par A(p) = T_s / T(p)
L'efficacité (ou efficiency) est définie par E(p) = A(p) / p.
A(p) doit être le plus près possible de p et E(p) le plus près possible de 1 (voire même plus grand).
Equation de Poisson
L'équation de Poisson en trois dimensions d'espace est notre premier problème modèle.Lien vers la page de l'exemple.
Equation de Helmholtz
L'équation de Helmholtz en trois dimensions d'espace est notre second problème modèle.Lien vers la page de l'exemple.
Equation de la Chaleur
L'équation de la Chaleur en trois dimensions d'espace est notre troisième problème modèle.Lien vers la page de l'exemple.
Comparaison des méthodes
L'équation de la Chaleur en trois dimensions d'espace sert de modèle pour comparer les performances des différentes méthodes proposées ci-dessus :- MPI décomposition de domaine avec MPI
- 1OMP partage du travail avec OpenMP (version 1 : parallélisation implicite)
- 2OMP partage du travail avec OpenMP (version 2 : parallélisation explicite)
- 3OMP partage du travail avec OpenMP (version 3 : décomposition de domaine)
- 1OMPI parallélisation à deux niveaux (version 1 : parallélisation implicite pour OpenMP)
- 2OMPI parallélisation à deux niveaux (version 2 : parallélisation explicite pour OpenMP)
Pour les différentes méthodes, on considère le même cas test et on étudie la scalabilité en augmentant le nombre de processus (ou threads) avec une masse de calcul identique par processus (ou threads).
Lien vers la page des comparaisons.
Conclusions et Tuning
Le choix de la technique de parallélisation doit aussi reposer sur l'architecture matérielle à disposition pour effectuer les calculs.La version MPI a l'avantage de tourner sur tout type de machine multi-processeur, qu'elle soit en une pièce (SGI Origin 2000, Compaq GS) ou sous forme de multi-noeuds (IBM SP, Compaq AlphaServer, cluster Linux, ...) alors que la version OpenMP nécessite une architecture à mémoire partagée. C'est le cas de la gamme Origin de SGI ou GS de Compaq et dans la limite de la taille des noeuds pour l'IBM SP (32 processeurs) et de l'AlphaServer de Compaq (4 processeurs), par exemple.
Les performances des codes dépendent d'une part de la parallélisation effectuée, du type des communications ou des synchronisations, du load-balancing, etc .... mais aussi de l'architecture sur laquelle les calculs sont effectués.
Une caractéristique importante est le débit (bande passante) entre la mémoire et les processeurs. Les noeuds IBM sont dotés de puces bi-processeurs qui se partagent un même cache L2 de 1.44 Mo. Deux processus qui tournent sur une même puce se partagent d'une part le cache L2 (soit 720 ko par processus) et d'autre part la bande passante entre ces deux processeurs et la mémoire. Ainsi ils ont chacun la moitié de la bande passante.
Si chaque processus est seul sur une puce, il dispose de l'intégralité de la bande passante soit le double du cas précédent.
Ainsi pour une application parallèle où les processus ont un comportement homogène, il peut y avoir concurrence pour accéder à la mémoire et donc détérioration des performances. Pour les exemples ci-dessus, les différences (constatées) peuvent aller jusqu'à 25%.
Des variables d'environnement que l'on peut paramétrer sont les suivantes (syntaxe propre au C-shell) :
Pour OpenMP :
setenv AIXTHREAD_SCOPE S
setenv SPINLOOPTIME 1000000
setenv YIELDLOOPTIME 1000000
AIXTHREAD_SCOPE indique que l'on souhaite mettre une seule thread par processeur physique.
SPINLOOPTIME indique au système que la thread doit rester active et verrouiller le processeur.
YIELDLOOPTIME complète l'action de SPINLOOPTIME.
Pour MPI :
setenv MP_SHARED_MEMORY yes
setenv MP_EUILIB ip
setenv MP_WAIT_MODE poll
MP_SHARED_MEMORY indispensable !! indique que les communications au sein des noeuds s'effectuent au travers de la mémoire (qui est partagée) et non pas au moyen de sockets ip.
MP_EUILIB indique que le job MPI ne peut pas être exécuté à cheval sur plusieurs noeuds (les performances du switch d'interconnection Colony ne sont pas à la hauteur). Les jobs resteront donc au sein d'un noeud.
MP_WAIT_MODE permet de spécifier le comportement d'un processus lorsqu'il est en attente (de réception par exemple) : le mode poll permet de faire avancer les processus ensemble, le mode yield est un mode d'attente active, on peut aussi mettre sleep dans le cas de processus maitre qui ne font strictement rien afin d'éviter de bloquer un processeur inutilement et de consommer du temps cpu (qui sera par conséquent décompté de l'attribution).