Résolution de grands systèmes linéaires

Techniques de renumérotation

   Une renumérotation des inconnues d’un système linéaire creux, c’est à dire une permutation de ses lignes et de ses colonnes, permet de limiter le phénomène de remplissage comme le montre l’exemple de la figure 1.2. Lors de la factorisation de la matrice de gauche, une structure pleine est créée. En renumérotant les équations en ordre inverse du plus grand au plus petit, on obtient la structure de droite pour laquelle il n’y a aucun remplissage au cours de la factorisation. Plusieurs techniques de renumérotation heuristiques peuvent être utilisées. C’est notamment le cas des variantes de l’algorithme de degré minimum (Tinney et Walker 1967; Liu 1985; Amestoy et al. 1996) ou celles de la dissection emboîtée (George 1973; George et Liu 1981). Le résultat obtenu est un arbre d’élimination représentant les dépendances dans les calculs et permettant de mener à bien les phases suivantes. Dans la pratique, cet arbre peut être construit directement à partir du graphe représentant la structure creuse du système.

Factorisation symbolique

   La factorisation symbolique de la matrice M consiste à déterminer le nombre et la position de ses nouveaux termes de remplissage qui seront créés au cours de la phase numérique. Cette procédure permet de construire la structure de données par blocs capable d’accueillir tous les coefficients de la matrice factorisée et donc d’optimiser le calcul effectif des coefficients des matrices triangulaires L et U. A priori, la factorisation symbolique de M requiert de faire toutes les itérations de l’algorithme d’élimination de Gauss en remplaçant le calcul des coefficients par un test sur leur pré-existence. Nous effectuons cela en incluant dans le super-arbre généré toutes les dépendances entre les inconnues du système linéaire.
Les super-arbres d’élimination et dépendances La construction d’un super-arbre d’élimination qui contient toutes les dépendances entre les inconnues n’est pas une chose aisée si nous partons des outils traditionnels de partitionnement de graphes. Par exemple, METIS fournit toute une variété de programmes qui peuvent être utilisés pour partitionner des graphes ou des maillages, ainsi que des programmes pour renuméroter les inconnues d’un système linéaire creux. Dans ces programmes, on utilise un algorithme de partitionnement à plusieurs niveaux en cherchant à chaque niveau, un séparateur de taille minimale qui découpe le graphe en deux sous-graphes de tailles égales. Mais aucun de ces programmes ne fournit des informations sur ces séparateurs dont nous avons besoin pour la construction et la recherche des dépendances entre les inconnues. En effet, les calculs effectués sur une inconnue i d’un séparateur peuvent dépendre de l’élimination d’une inconnue j d’un supernœud situé à un niveau inférieur, si le terme mij est non nul. Nous voulons donc trouver pour chacun des super-nœuds dans l’arbre, toutes les inconnues qui sont en relation avec les siennes et celles de ses descendants. Nous qualifierons les inconnues propres à un super-nœud d’inconnues internes et celles qui sont dépendantes de ce super-nœud seront appelées inconnues interfaces. Pour construire les super-arbres d’élimination renfermant toutes les dépendances, nous proposons donc deux approches qui dépendent du type de découpage choisi. Dans la première approche, le découpage est effectué sur les éléments finis d’un maillage et il est géré par le partitionneur METIS. Cela ne nous permet pas de conserver toute l’historique sur la partition des éléments. Le nombre de sous-structures obtenues dans cette partition peut être une puissance de deux ou bien être quelconque. Dans la deuxième approche, nous effectuons directement le découpage sur les inconnues du système, et non plus sur les éléments du maillage. Cela assure le respect de l’ordre initial des inconnues, donc de la structure de la matrice. En découpant, nous intervenons dans les sous-programmes de METIS pour récupérer à chaque niveau les séparateurs trouvés. Le nombre de sous-structures présentes dans ce deuxième type de découpage ne s’écrit que sous la forme d’une puissance de deux. La première approche consiste, en même temps, à reconstituer toute l’historique sur la partition du maillage éléments finis composée d’un nombre ndoms de sous-structuress Ii et à construire le super-arbre d’élimination contenant toutes les dépendances existantes entre les inconnues. Pour cela, nous mettons d’abord en place, grâce à l’algorithme 1.1, la structure de l’arbre en calculant sa hauteur h et le nombre nbnode(l) de super-nœuds situés à chacun de ses niveaux l. Ensuite, nous cherchons les nœuds (ou inconnues) internes et les nœuds (ou inconnues) interfaces de chacun des sous-domaines Ii que nous placerons au plus bas niveau de  l’arbre (l = 0)

Évaluation du nombre optimal de sous-structures

   Afin de déterminer le nombre optimal de sous-structures engendrées par bissection récursive et conduisant à de meilleures performances du solveur, nous résolvons d’abord des problèmes de tailles neqs différentes (41472, 107811 et 397953 équations), où chacun d’eux est partitionné en un nombre ndoms de sous-domaines (64, 128, 256, 512 et 1024). Ensuite, nous analysons les temps de calcul nécessaires pour leur résolution. Mais auparavant, nous porterons une attention particulière sur la qualité des découpages et leurs effets sur les temps écoulés dans la phase numérique. Les tableaux 1.2, 1.3 et 1.4 présentent respectivement les résultats sur les découpages des problèmes de tailles 41472, 107811 et 397953. D’une part, ces résultats portent sur le nombre total sepsz d’inconnues internes aux séparateurs, le pourcentage sepsz/neqs, la taille mid (total des inconnues internes et interfaces) des sous-structures les plus représentées, la taille moyenne avrg de toutes les sous-structures, et le rapport ratio entre mid et avrg. D’autre part, ils nous informent sur les temps nécessaires pour réaliser la phase numérique : le temps tloc de factorisation des matrices locales associées aux sous-structures et celui tschur de construction et de factorisation du complément de Schur.

Solveurs destinés aux architectures à mémoire distribuée

   L’évolution des architectures des calculateurs scientifiques dans les dernières années a consacré l’usage systématique du parallélisme multiprocesseurs à mémoire distribuée pour les grands systèmes linéaires. Mémoire distribuée veut dire que chaque processeur (ou ensemble de processeurs) dispose d’une mémoire locale et ne peut accéder directement qu’à celle-ci. Les données résidant dans les différentes mémoires locales ne peuvent s’échanger que par des transferts réalisés par un réseau d’interconnexion (figure 2.1). Le parallélisme à mémoire distribuée offre une solution de conception pour certains solveurs directs creux parmi lesquels nous pouvons citer les codes PSPASES (Gupta et al. 1994) et DSCPack (Raghavan 2001) qui sont spécialisés dans les systèmes symétriques définis positifs (SDP), et SuperLU_Dist (Li et Demmel 2003) pour des systèmes non symétriques. Présentement, la seule bibliothèque parallèle adaptée à ces deux types de systèmes est MUMPS (Amestoy et al. 2000). PSPASES est l’un des premiers solveurs parallèles développé pour les systèmes SDP qui a réussi à résoudre en un temps réduit un problème comportant un million d’équations sur 256 processeurs d’un CRAY T3E. L’algorithme parallèle dans SuperLU_DIST est basé sur une partition en super-nœuds renfermant un parallélisme de haut niveau dans la mise à jour du complément de Schur. Des algorithmes d’élimination sur des graphes orientés acycliques sont utilisés pour identifier la dépendance entre les super-nœuds. Le modèle de programmation parallèle utilisé dans ces solveurs est basé sur les librairies d’échanges de messages (MPI), ce qui les rend donc portables sur les calculateurs parallèles à mémoire distribuée.

Le concept général

   La notion de parallélisme à mémoire partagée, fréquemment appelé multithreading, est rendue nécessaire par l’utilisation qui est faite aujourd’hui des processeurs multi-cœurs. Dans cette approche, chaque thread a un accès partagé à la mémoire globale du processeur et les échanges de données entre les différents threads se font simplement par des lectures/écritures vers la même zone mémoire cache. Il faut s’assurer qu’aucune tâche ne peut lire une donnée qui n’est que partiellement mise à jour du fait que les opérations peuvent utiliser plusieurs cycles d’horloge. Si une partie de la cohérence des données peut être assurée par le matériel, la totalité de la consistance des données reste à la charge du développeur. La parallélisation des applications n’étant pas toujours une chose aisée, des outils et des compilateurs spécialisés ont vu le jour pour faciliter le multithreading. Dans ce cas, le parallélisme est géré directement par ces outils. Le compilateur génère alors un code parallèle, utilisant le multi threading en suivant les indications données par le programmeur.

Le standard OpenMP

   OpenMP est une interface standard de haut niveau pour une programmation parallèle sur machine à mémoire partagée permettant une gestion de la distribution des calculs entre plusieurs processeurs multi-cœurs. Basée sur les techniques du multi-threading, OpenMP peut être considérée comme l’un des grands standards au service du calcul scientifique. Cette interface de programmation est supportée sur la plupart des plate-formes, incluant Unix et Windows, pour les langages de développement Fortran, C et C++. Grâce à des directives insérées par le développeur dans le code et grâce à une analyse des dépendances entre les données, un compilateur OpenMP est capable d’extraire du parallélisme des boucles de calcul. Les directives ne sont que des commentaires activables par le compilateur grâce une option adéquate. Ainsi, un compilateur ne supportant pas OpenMP les ignorera. C’est ce qui permet de garantir la portabilité des codes OpenMP et la maintenance d’une version unique, séquentielle et parallèle du code. La parallélisation des portions de code indépendantes va être construite par le compilateur en utilisant des threads. Ce standard est donc une approche simplifiée pour paralléliser des codes via des threads. Certaines librairies mathématiques sont déjà optimisées avec OpenMP pour s’exécuter sur tous les cœurs d’un nœud de calcul. C’est le cas, par exemple, de la librairie MKL d’Intel qui regroupe un ensemble de sous-programmes scientifiques telles que BLAS, LAPACK. Des variables d’environnement permettent de contrôler l’environnement d’exécution comme par exemple le nombre de threads. Par défaut, le nombre de threads utilisés dans BLAS est égal au nombre de cœurs de la machine. Dans certains cas, pour être optimal, il faut ajuster le nombre de threads à l’exécution avec la variable d’environnement OMP_NUM_THREADS ou bien faire appel à la fonction omp_set_num_threads().

Les temps de calcul pour la phase d’analyse

   Les temps nécessaires pour analyser le système linéaire à résoudre concernent la renumérotation des inconnues (bissection récursive et construction du super-arbre d’élimination) et la factorisation symbolique des associées aux sous-structures. Sur la figure 2.15, nous présentons l’accélération des calculs quand le nombre de threads créés croît. Nous remarquons que l’accélération de la version optimisée avec OpenMP est, sans surprise, identique à 1 quel que soit le nombre de threads. Les temps obtenus en parallèle sont identiques à celui obtenu en séquentiel. La raison est que nous n’utilisons pas des sousprogrammes de BLAS dans cette phase d’analyse. Nous observons aussi une faible évolution de l’accélération de la version Pthreads. Elle passe de 1 à 1.07 quand le nombre de threads passe de 1 à 8. Cette faible amélioration des temps de calcul est due au fait que seule l’étape de factorisation symbolique (SymbFact) était parallélisable. Or, cette étape ne coûte pas grand chose comparée à celle de la construction de l’arbre (BuildTree) comme nous pouvons l’observer sur le tableau 2.1. Les bonnes performances en parallèle obtenues pendant la factorisation symbolique (figure 2.16) sont donc fortement influencées par l’étape de bissection récursive sur les inconnues.

Table des matières

Liste des figures
Liste des tableaux
Liste des algorithmes
Introduction générale
1 Mise en œuvre d’un solveur direct 
1.1 Le cadre général
1.1.1 Matrices creuses et remplissage
1.1.2 Techniques de renumérotation
1.1.3 Arbre d’élimination
1.2 La phase d’analyse 
1.2.1 Dissection emboîtée
1.2.2 Factorisation symbolique
1.3 La phase numérique
1.3.1 Mise à l’échelle ou « scaling » de la matrice originale
1.3.2 Factorisation numérique
1.3.3 Condensation statique
1.4 La résolution des systèmes triangulaires 
1.4.1 Phase de descente
1.4.2 Résolution du problème condensé
1.4.3 Phase de remontée
1.5 La prise en compte des systèmes non inversibles 
1.5.1 Détection locale des singularités
1.5.2 Traitement des modes à énergie nulle
1.5.3 Illustration du déroulement de l’algorithme
1.6 L’évaluation des performances 
1.6.1 Comparaison des approches proposées dans la phase d’analyse
1.6.2 Évaluation du nombre optimal de sous-structures
1.6.3 Comparaison des temps d’exécution
1.6.4 Comparaison des temps de descente-remontée
1.6.5 Analyse de la complexité du solveur
Conclusion
2 Parallélisation du solveur 
2.1 L’état de l’art
2.1.1 Solveurs destinés aux architectures à mémoire distribuée
2.1.2 Solveurs adaptés aux architectures à mémoire partagée
2.1.3 Solveurs directs basés sur une programmation hybride
2.1.4 Discussions
2.2 La programmation multi-cœurs
2.2.1 Motivation
2.2.2 Architectures multi-cœurs
2.2.3 Parallélisme à mémoire partagée : le multi-threading
2.3 Les démarches pour la parallélisation du solveur
2.3.1 Version multi-threads avec Pthreads
2.3.2 Version multi-threads avec OpenMP
2.3.3 Évaluation des deux versions multi-threads
2.3.4 Mise en place d’une version multi-threads hybride
2.4 Les performances en multi-threading
2.4.1 Analyse des performances
2.4.2 Comparaison avec DSCPack et MUMPS
Conclusion
3 Amélioration et évaluation des performances de la méthode FETI 
Introduction
3.1 La présentation de la méthode FETI 
3.1.1 Description de la version initiale de FETI
3.1.2 Résolution itérative du problème d’interface
3.1.3 Algorithme détaillé de FETI-1
3.2 Les facteurs critiques d’une méthode FETI
3.2.1 Choix d’un solveur local
3.2.2 Traitement des singularités et calcul des pseudo-inverses
3.2.3 Gestion du problème grossier
3.3 L’intégration d’un parallélisme à deux niveaux dans FETI
3.3.1 Extraction en parallèle des colonnes de B (s)
3.3.2 Calcul en parallèle des pseudo-inverses et pseudo-solutions
3.4 L’évaluation des performances de la méthode FETI
3.4.1 Environnement d’évaluation
3.4.2 Critères de mesure de performances
3.4.3 Description des cas tests parallèles
3.4.4 Analyse de l’efficacité de FETI
3.4.5 Analyse de l’extensibilité de FETI
Conclusion
Conclusion générale
A Annexes
A.1 Factorisation par élimination de Gauss
A.1.1 Les généralités de la méthode
A.1.2 Les inconvénients de la méthode
A.2 Algorithmes de renumérotation
A.2.1 L’algorithme de dissection emboîtée
A.2.2 L’algorithme de degré minimum et ses variantes
A.2.3 L’algorithme de Cuthill – McKee
A.3 Présentation de la méthode des éléments finis
A.3.1 Le problème modèle
A.3.2 La mise sous forme variationnelle d’un problème d’EDP
A.3.3 Le choix d’un maillage et discrétisation
A.3.4 Le problème sous forme matricielle
A.4 Analyse des performances des versions multi-threads du solveur Dissection
A.4.1 Les temps de calcul pour la phase d’analyse
A.4.2 Les temps de calcul pour la phase numérique
A.4.3 Les temps de calcul pour la phase de descente-remontée
A.4.4 Les temps d’exécution en fonction du nombre de threads
A.4.5 Les conclusions
Bibliographie
Notations

Télécharger le rapport complet