9 · Solveurs et algèbre linéaire¶
Le contenu mathématique que le cursus FISA ne donne pas, parce qu'AEDP24 est inaccessible et qu'ANAF23 appartient au parcours Mathématiques appliquées. C'est aussi la matière des deux benchmarks qui classent les supercalculateurs.
Difficulté : ★★★★ · ⏱ 40 h
Chapitre optionnel, à décider consciemment
C'est le chapitre le plus long et le plus lent à rentabiliser du socle. Si votre objectif est une compétition dans dix-huit mois, passez-le : lisez uniquement les sections « Les niveaux de BLAS » et « HPL et HPCG », qui sont directement opérationnelles, et gardez le reste pour un M2 ou une thèse.
Si votre objectif est une carrière en calcul scientifique — CEA, EDF, ONERA, météorologie, recherche — ce chapitre est au contraire le plus important du socle, parce que les gains les plus considérables sont algorithmiques et non techniques.
Pourquoi c'est central¶
Une part écrasante du temps de calcul mondial se passe à résoudre des systèmes linéaires, directement ou indirectement. Toute discrétisation d'équation aux dérivées partielles — donc toute simulation de fluide, de structure, de champ électromagnétique, de transfert thermique — se ramène à cela.
Et c'est aussi, littéralement, ce que mesurent les classements :
- HPL (High Performance Linpack), qui définit le Top500, résout un système linéaire dense par factorisation LU avec pivotage partiel ;
- HPCG (High Performance Conjugate Gradient), le second classement officiel, résout un système creux par gradient conjugué préconditionné multigrille.
Les deux donnent des classements très différents, et comprendre pourquoi est le meilleur résumé possible de tout ce chapitre.
Les niveaux de BLAS, et pourquoi cela structure tout¶
Le point le plus opérationnel du chapitre. Les Basic Linear Algebra Subprograms sont classés en trois niveaux :
| Niveau | Opération type | Données | Opérations | Intensité |
|---|---|---|---|---|
| 1 | \(y \leftarrow \alpha x + y\) (AXPY), produit scalaire | \(O(n)\) | \(O(n)\) | \(O(1)\) — faible |
| 2 | \(y \leftarrow \alpha A x + \beta y\) (GEMV) | \(O(n^2)\) | \(O(n^2)\) | \(O(1)\) — faible |
| 3 | \(C \leftarrow \alpha A B + \beta C\) (GEMM) | \(O(n^2)\) | \(O(n^3)\) | \(O(n)\) — croissante |
La conséquence est immense. Les niveaux 1 et 2 ont une intensité arithmétique constante et faible : ils sont irrémédiablement limités par la bande passante mémoire, et leur performance plafonne à quelques pour cent du pic de calcul, quelle que soit la qualité de l'implémentation.
Le niveau 3 a une intensité qui croît avec la taille : on peut le bloquer pour le cache et atteindre 85 à 95 % du pic.
D'où le principe directeur de toute l'algorithmique numérique haute performance : reformuler les algorithmes pour exprimer le maximum du travail en BLAS 3.
C'est exactement ce que fait LAPACK par rapport à son prédécesseur LINPACK : les mêmes algorithmes mathématiques, réécrits par blocs pour que le travail dominant soit du GEMM. Le gain n'est pas de quelques pour cent, c'est un ordre de grandeur.
Le réflexe à acquérir
Devant un algorithme numérique, la question à se poser est : quelle fraction du travail est en BLAS 3 ? Si la réponse est « aucune », c'est un algorithme limité par la mémoire et sa performance sera médiocre par nature. Si la réponse est « l'essentiel », il peut approcher le pic.
Corollaire pratique : n'écrivez jamais votre propre GEMM en production. OpenBLAS, BLIS et MKL y ont consacré des années-homme. Appelez-les. Écrire un GEMM soi-même est un exercice pédagogique (et excellent), pas une décision d'ingénierie.
Les méthodes directes¶
On factorise la matrice, puis on résout par substitutions.
| Méthode | Condition sur \(A\) | Coût | Usage |
|---|---|---|---|
| LU avec pivotage partiel | Générale | \(\frac{2}{3}n^3\) | Le cas général. C'est HPL |
| Cholesky (\(A = LL^\top\)) | Symétrique définie positive | \(\frac{1}{3}n^3\) | Deux fois moins cher, pas de pivotage |
| QR | Générale | \(\frac{4}{3}n^3\) | Moindres carrés, plus stable |
| LDL\(^\top\) | Symétrique | \(\frac{1}{3}n^3\) | Symétrique non définie positive |
La formule à connaître : le coût de LU est \(\frac{2}{3}n^3\) opérations flottantes. C'est elle qui permet de convertir un temps HPL en FLOPS :
où \(n\) est l'ordre de la matrice et \(t\) le temps en secondes. Le terme \(2n^2\) vient des substitutions. C'est la formule exacte utilisée par HPL.
En creux, les méthodes directes souffrent du remplissage (fill-in) : la factorisation d'une matrice creuse produit des facteurs beaucoup moins creux que la matrice initiale. Le coût mémoire peut devenir prohibitif. D'où l'importance du réordonnancement (algorithmes de dissection emboîtée, nested dissection, et de degré minimum) qui minimise le remplissage — c'est un problème de théorie des graphes, et c'est ce que font METIS et SCOTCH (ce dernier étant développé à Bordeaux).
Bibliothèques : MUMPS (français, Inria et CERFACS), SuperLU, PARDISO, UMFPACK.
Les méthodes itératives¶
On construit une suite qui converge vers la solution, sans jamais factoriser. Chaque itération coûte essentiellement un produit matrice-vecteur.
Les méthodes de Krylov, par ordre de généralité :
| Méthode | Condition sur \(A\) | Remarque |
|---|---|---|
| Gradient conjugué (CG) | Symétrique définie positive | Optimal dans ce cas. C'est HPCG |
| MINRES | Symétrique | — |
| BiCGStab | Générale | Peu de mémoire, convergence parfois erratique |
| GMRES | Générale | Robuste, mais mémoire croissante — d'où le redémarrage GMRES(\(m\)) |
Le résultat théorique à connaître pour le gradient conjugué : le nombre d'itérations nécessaires pour atteindre une précision donnée croît comme la racine carrée du conditionnement :
où \(\kappa(A)\) est le conditionnement de \(A\) (rapport de la plus grande à la plus petite valeur propre) et \(\varepsilon\) la précision relative visée.
C'est pourquoi le préconditionnement est décisif. On résout \(M^{-1}Ax = M^{-1}b\) avec \(M\) choisi de sorte que \(\kappa(M^{-1}A) \ll \kappa(A)\). Un bon préconditionneur peut faire passer le nombre d'itérations de dix mille à cinquante.
Les préconditionneurs, du plus simple au plus puissant :
| Préconditionneur | Coût | Efficacité | Parallélisable ? |
|---|---|---|---|
| Jacobi (diagonal) | Négligeable | Faible | Parfaitement |
| Gauss-Seidel, SOR | Faible | Moyenne | Mal (séquentiel par nature) |
| Factorisation incomplète (ILU) | Moyen | Bonne | Mal |
| Schwarz additif par blocs | Moyen | Bonne | Bien |
| Multigrille (géométrique ou algébrique) | Élevé | Excellente | Bien, avec du soin |
Le multigrille est la famille de méthodes la plus importante du domaine : pour les problèmes elliptiques, elle atteint une complexité optimale, c'est-à-dire un coût linéaire en nombre d'inconnues, avec un nombre d'itérations indépendant de la taille du problème. C'est un résultat remarquable, et c'est ce qui permet de résoudre des systèmes à des milliards d'inconnues.
Le préconditionneur de HPCG est précisément un multigrille symétrique par Gauss-Seidel.
Bibliothèques : hypre (Lawrence Livermore, avec le multigrille algébrique BoomerAMG), PETSc, Trilinos, AMGX (NVIDIA, sur GPU).
Direct contre itératif : l'arbitrage¶
La règle de décision
| Critère | Direct | Itératif |
|---|---|---|
| Matrice dense | Oui | Non |
| Matrice creuse, 2D | Souvent bon | Bon |
| Matrice creuse, 3D | Remplissage prohibitif | Oui |
| Plusieurs seconds membres | Oui (une factorisation, N résolutions) | Coûteux |
| Mal conditionnée | Oui (robuste) | Difficile sans bon préconditionneur |
| Très grande taille | Mémoire prohibitive | Oui |
| Prévisibilité du temps | Oui (coût connu à l'avance) | Non (dépend de la convergence) |
En pratique, les gros codes 3D utilisent presque tous des méthodes itératives préconditionnées, et souvent un solveur direct sur les sous-problèmes locaux (comme préconditionneur de type Schwarz).
Pourquoi HPL et HPCG donnent des classements différents¶
La question la plus instructive du chapitre.
| HPL | HPCG | |
|---|---|---|
| Problème | Système dense, factorisation LU | Système creux, CG préconditionné multigrille |
| Travail dominant | GEMM, donc BLAS 3 | SpMV, donc BLAS 2 creux |
| Intensité arithmétique | Élevée et croissante | Très faible |
| Régime | Compute bound | Memory bound |
| Communications | Volumineuses mais recouvrables | Allreduce à chaque itération, synchronisant |
| Performance atteinte | 60 à 90 % du pic | 1 à 5 % du pic |
| Représentativité des codes réels | Faible | Élevée |
La conclusion : HPL mesure la capacité de calcul brute, HPCG mesure ce qu'une application réelle obtient. L'écart entre les deux, sur la même machine, est typiquement d'un facteur vingt à cinquante. C'est la mesure de l'écart entre les machines telles qu'elles sont conçues et les codes tels qu'ils sont.
HPCG a été créé précisément pour cela, par Dongarra, Heroux et Luszczek : pour donner aux constructeurs une incitation à améliorer la bande passante et la latence des communications, et pas seulement le pic de calcul.
Les bibliothèques à connaître¶
| Bibliothèque | Nature | À savoir |
|---|---|---|
| BLAS / LAPACK | Interfaces de référence | Les API, pas les implémentations. Tout en dépend |
| OpenBLAS, BLIS | Implémentations libres de BLAS | BLIS a une architecture plus claire et documentée |
| Intel MKL, AMD AOCL, Arm PL | Implémentations vendeur | Souvent les plus rapides sur leur matériel |
| ScaLAPACK | LAPACK distribué (MPI) | Ancien, encore utilisé. Grille de processus 2D |
| SLATE | Successeur moderne de ScaLAPACK | Conçu pour les architectures hétérogènes |
| PETSc | Cadre complet de solveurs | Le plus important à connaître : solveurs, préconditionneurs, maillages, TS, optimisation |
| Trilinos | Équivalent, plus modulaire, C++ | Alternative américaine à PETSc |
| hypre | Préconditionneurs, dont BoomerAMG | La référence en multigrille algébrique |
| MUMPS | Solveur direct creux distribué | Français, de très bonne qualité |
| SuiteSparse | Solveurs creux séquentiels (UMFPACK, CHOLMOD) | La référence en séquentiel |
| FFTW | Transformée de Fourier | L'implémentation de référence ; à connaître pour son mécanisme de « plans » |
| METIS, SCOTCH, ParMETIS | Partitionnement de graphes | Pour le réordonnancement et la décomposition de domaine. SCOTCH est bordelais |
| cuBLAS, cuSPARSE, cuSOLVER, AMGX | Équivalents GPU | — |
PETSc, l'investissement le plus rentable du chapitre
Si vous ne devez apprendre qu'une bibliothèque, apprenez PETSc. Raisons :
- elle couvre tout : vecteurs et matrices distribués, une trentaine de solveurs de Krylov, une cinquantaine de préconditionneurs, des maillages structurés et non structurés, l'intégration temporelle, l'optimisation ;
- elle permet de changer de solveur et de préconditionneur par options de
ligne de commande, sans recompiler.
-ksp_type gmres -pc_type hypre. Cela transforme l'expérimentation numérique ; - elle est utilisée par un très grand nombre de codes de production ;
- sa documentation et ses tutoriels sont exemplaires, et son équipe est réactive.
Le coût d'entrée est de vingt à trente heures. Le retour est de pouvoir résoudre efficacement n'importe quel système linéaire ou non linéaire creux pour le reste de votre carrière.
Exercices¶
E1 · ★★ ⏱ 3 h — Les trois niveaux de BLAS. Mesurer les GFLOPS de daxpy,
dgemv et dgemm d'OpenBLAS, en balayant la taille. Tracer les trois courbes sur
le même graphique avec le pic de la machine en pointillés. Constater que seul le
GEMM s'approche du pic, et vérifier que le rapport correspond aux intensités
arithmétiques calculées.
E2 · ★★★ ⏱ 5 h — LU à la main, et par blocs. Implémenter une factorisation LU
avec pivotage partiel en trois versions : non bloquée (BLAS 1 et 2), bloquée avec
appels à dgemm pour les mises à jour, et par appel direct à dgetrf de LAPACK.
Mesurer les trois. L'écart entre la première et la deuxième démontre le principe du
BLAS 3 ; l'écart avec LAPACK mesure ce que vaut une décennie d'optimisation.
E3 · ★★★ ⏱ 4 h — Gradient conjugué et conditionnement. Implémenter un gradient conjugué, et le tester sur une famille de matrices de conditionnement croissant. Tracer le nombre d'itérations en fonction de \(\sqrt{\kappa}\) et vérifier la relation théorique. Puis ajouter un préconditionneur de Jacobi et retracer.
E4 · ★★★ ⏱ 4 h — SpMV, trois formats. Implémenter un produit matrice-vecteur creux dans trois formats de stockage : COO, CSR, et ELLPACK ou un format par blocs. Mesurer sur plusieurs matrices de structures différentes (issues de la collection SuiteSparse, qui est publique et standard). Constater que le meilleur format dépend de la matrice. Puis comparer à cuSPARSE si vous avez un GPU.
E5 · ★★★ ⏱ 5 h — PETSc, premier contact. Résoudre un problème de Poisson 3D avec PETSc, puis balayer les solveurs et préconditionneurs par options de ligne de commande : CG, GMRES, BiCGStab, avec Jacobi, ILU, Schwarz additif et hypre. Tracer le temps total et le nombre d'itérations pour chaque combinaison. Vous obtiendrez une matrice de résultats instructive, et vous comprendrez pourquoi PETSc est conçu ainsi.
E6 · ★★★ ⏱ 4 h — HPL, pour de vrai. Installer et régler HPL sur une machine ou
un petit cluster. Les paramètres du fichier HPL.dat à comprendre : N (ordre de
la matrice, à choisir pour occuper 80 à 90 % de la mémoire disponible), NB
(taille de bloc, typiquement entre 100 et 400, à ajuster), P et Q (grille de
processus, avec \(P \times Q\) = nombre de rangs et \(P \le Q\) de préférence).
Balayer les paramètres et trouver l'optimum. Comparer au pic théorique.
Cet exercice est directement une épreuve de compétition, et il prend plus de temps qu'on ne croit : c'est précisément pour cela qu'il faut l'avoir fait avant.
E7 · ★★★★ ⏱ 6 h — HPCG, et l'écart. Installer et exécuter HPCG sur la même
machine que E6. Comparer les deux chiffres et expliquer le facteur d'écart par
l'analyse d'intensité arithmétique et par le coût de l'Allreduce. Rédiger deux
pages. C'est un excellent sujet d'exposé, et un excellent sujet pour l'évaluation
par lecture d'article de PDSP35.
E8 · ★★★★ ⏱ 6 h — Multigrille à la main. Implémenter un multigrille géométrique à deux puis à plusieurs niveaux pour un problème de Poisson 1D ou 2D : lissage par Jacobi ou Gauss-Seidel, restriction, prolongation, cycle en V. Mesurer le nombre d'itérations en fonction de la taille du problème et constater qu'il ne croît pas. C'est l'un des plus beaux résultats du calcul scientifique, et le constater soi-même est marquant.
Ressources¶
Priorité 1 :
- Lloyd Trefethen, David Bau, Numerical Linear Algebra, SIAM, 1997. Quarante courtes leçons, remarquablement écrites. Le meilleur point d'entrée, et de loin.
- Yousef Saad, Iterative Methods for Sparse Linear Systems, 2ᵉ éd., SIAM, 2003. PDF gratuit sur le site de l'auteur. La référence sur les méthodes de Krylov et le préconditionnement.
- La documentation et les tutoriels de PETSc (
petsc.org). Exemplaires.
Priorité 2 :
- Gene Golub, Charles Van Loan, Matrix Computations, 4ᵉ éd., Johns Hopkins, 2013. La référence encyclopédique. À consulter.
- James Demmel, Applied Numerical Linear Algebra, SIAM, 1997. Plus orienté algorithmique et stabilité.
- Barrett et al., Templates for the Solution of Linear Systems, SIAM, gratuit : un recueil compact des algorithmes de Krylov en pseudo-code. Très pratique pour implémenter.
- Briggs, Henson, McCormick, A Multigrid Tutorial, 2ᵉ éd., SIAM. Court et clair sur le multigrille.
- Le cours CS267 de Berkeley (Demmel), dont plusieurs leçons portent exactement sur ce chapitre.
Priorité 3 :
- Goto & van de Geijn, « Anatomy of High-Performance Matrix Multiplication », ACM TOMS, 2008, et les articles de BLIS (van Zee & van de Geijn, 2015).
- Les spécifications de HPL (
netlib.org/benchmark/hpl/) et de HPCG (hpcg-benchmark.org), y compris l'article de Dongarra, Heroux et Luszczek sur la motivation de HPCG. - La collection SuiteSparse Matrix Collection (anciennement University of Florida) : des milliers de matrices creuses réelles, standard pour les comparaisons.
- Sur les méthodes communication-avoiding : les travaux de Demmel, Ballard, Hoemmen, Grigori. Le plus élégant du domaine.
- MUMPS (
mumps-solver.org) et SCOTCH, deux projets français de référence.
À retenir¶
Les cinq idées du chapitre
- Les niveaux de BLAS structurent tout. Le niveau 3 (GEMM) peut approcher le pic ; les niveaux 1 et 2 sont limités par la mémoire. L'algorithmique numérique haute performance consiste à reformuler en BLAS 3.
- N'écrivez jamais votre propre GEMM en production. Appelez OpenBLAS, BLIS ou MKL.
- Direct pour le dense et les seconds membres multiples, itératif préconditionné pour le creux 3D de grande taille.
- Le préconditionnement décide de tout en itératif, et le multigrille est optimal pour les problèmes elliptiques : nombre d'itérations indépendant de la taille.
- L'écart HPL / HPCG, d'un facteur vingt à cinquante, mesure l'écart entre la machine et les applications réelles. Savoir l'expliquer est le meilleur résumé possible de la compréhension du domaine.
Et la formule opérationnelle
Le coût d'une factorisation LU est \(\frac{2}{3}n^3\) opérations flottantes. C'est elle qui convertit un temps HPL en FLOPS, et vous l'utiliserez souvent.
Chapitre suivant : Énergie et efficacité.