# Optimiser une simulation SPH, partie 1 : la vitesse offerte par le compilateur Avant de toucher la moindre ligne de ma simulation de fluide, j'ai laissé le compilateur faire le travail. Changer simplement les flags m'a fait passer de 1x à environ 40x. ## La base de référence J'ai écrit une simulation de fluide SPH simple basée sur « Particle-Based Fluid Simulation for Interactive Applications » de Matthias Müller, David Charypar et Markus Gross. ![La simulation SPH non optimisée, tournant péniblement lentement](./baseline-sim.png) Le projet n'est pas optimisé par défaut. C++20, `std::vector` partout, `glm` pour les maths, précision `double`, build de debug. Ça tourne. C'est juste douloureusement lent. Je me suis alors demandé combien de performance je pouvais tirer de ma machine (un ThinkPad P14s Gen 5 avec un AMD Ryzen 8 PRO 8840HS) pour cette simulation. Dans cette série, je vais parcourir tout le cheminement. Cette première partie porte sur le gain le moins coûteux disponible : laisser le compilateur optimiser mon code naïf avant d'en réécrire quoi que ce soit. ## Le flag -O > Des flags contrôlant le niveau d'optimisation à effectuer. Derrière ce petit flag se cache une énorme quantité de machinerie clang. Il finit par être transmis au [`PassBuilder`](https://llvm.org/doxygen/classllvm_1_1PassBuilder.html) de LLVM, et si vous jetez un œil à `llvm::OptimizationLevel` vous obtenez un résumé assez honnête de chaque niveau : - **0** : désactiver autant d'optimisations que possible. - **1** : optimiser rapidement sans détruire la capacité de débogage. - **2** : optimiser pour une exécution rapide sans faire exploser le temps de compilation ou la taille du code. - **3** : optimiser pour une exécution rapide, un point c'est tout. Ça, c'est la théorie. Voyons ce qui sort vraiment à l'autre bout. ## Lire l'assembleur : O0 vs O1 À `-O0` le code généré pour ma fonction `compute_forces` est exactement ce à quoi on s'attendrait : chaque petite chose est écrite explicitement, rien n'est réutilisé, et les opérateurs `glm` sont de vrais appels de fonction. Voici un extrait représentatif, calculant la distance entre deux particules : ```asm lea rax, [rbp-0x40] mov rsi, rcx mov rdi, rax call glm::operator- lea rax, [rbp-0x40] mov rdi, rax call glm::length<3, double, 0> movq rax, xmm0 mov qword [rbp-0x178 {r}], rax movsd xmm0 {h}, qword [rbp-0x1c0] comisd xmm0, qword [rbp-0x178 {r}] ``` Deux appels de fonction juste pour soustraire deux vecteurs et calculer une longueur. Voici maintenant la même logique à `-O1` : ```asm movsd xmm1, qword [rbx+0x10] subsd xmm1, qword [r12+0x10] movsd xmm0, qword [rbx+0x8] subsd xmm0, qword [r12+0x8] movsd xmm2, qword [rbx] subsd xmm2, qword [r12] mulsd xmm0, xmm0 mulsd xmm2, xmm2 addsd xmm0, xmm2 mulsd xmm1, xmm1 addsd xmm0, xmm1 sqrtsd xmm0, xmm0 movsd xmm3 {smoothing_radius}, qword [rsp] comisd xmm3, xmm0 ``` Quelques choses notables se sont produites. ### Inlining Les appels à `glm::operator-` et `glm::length` ont disparu. Le compilateur a tiré leurs corps directement dans la boucle, donc il ne reste plus que des soustractions, des multiplications, une addition, et un `sqrtsd` Un `call` n'est pas gratuit. Il faut préparer les arguments, sauter quelque part, puis revenir. Cela impacte le cache d'instructions et le pipeline. L'inlining élimine tout ça. Le compromis, c'est que le binaire devient plus gros, puisque le corps est désormais dupliqué partout où il était appelé. ### Range-for abaissé en arithmétique de pointeurs À `-O0`, ma boucle `for (auto& p : particles)` est prise très littéralement. Le compilateur matérialise `__for_range`, `__for_begin` et `__for_end` en véritables emplacements de pile, et appelle `begin()` et `end()` : ```asm call std::vector::begin mov qword [rbp-0x1e0 {__for_begin}], rax call std::vector::end mov qword [rbp-0x1d8 {__for_end}], rax ``` À `-O1` tout ça disparaît. Il se contente de charger `_M_start` et `_M_finish` directement depuis le vecteur, fait avancer le pointeur de particule courant dans un registre, et l'incrémente de la taille d'une `Particle` à chaque itération : ```asm add r12, 0x58 cmp qword [rsp+0x78], r12 ``` On se rapproche de la boucle classique sur tableau en C. L'abstraction est toujours là dans le code source, elle ne coûte simplement rien à l'exécution. ### Pas d'autovectorisation Une chose qui ne s'est *pas* produite : la vectorisation. L'auto-vectorisation, c'est le compilateur qui prend une boucle que vous avez écrite comme si elle traitait un élément à la fois, et la réécrit pour en traiter plusieurs à la fois grâce aux instructions SIMD du CPU. Un simple `addsd` additionne un `double`. La version vectorielle `addpd` en additionne deux (ou quatre, ou huit) empaquetés côte à côte en une seule instruction. Même coût, plus de travail effectué. Les boucles ici sont restées scalaires. Chaque opération est sous sa forme à élément unique. Je reviendrai sur la raison dans une partie ultérieure, car il s'avère que la disposition des données se bat ici contre le compilateur. ## O2 et O3 Ça change pas mal l'assembleur, et j'avoue que je ne peux pas en lire chaque ligne. Mais voici ce que j'ai pu en tirer. À `-O2`, les accumulateurs vivent dans des registres pendant toute la boucle interne. À `-O1` les six composantes de force (`pressure_force.x/y/z` et `viscous_force.x/y/z`) résidaient dans des emplacements de pile fixes et étaient rechargées et stockées à chaque itération. À `-O2` elles sont maintenues dans des registres tout au long de la boucle et ne touchent la mémoire que quand c'est vraiment nécessaire. La boucle interne devient donc principalement des `mulsd`/`divsd`/`addsd` registre-à-registre, avec très peu de trafic sur la pile. L'autre chose : le chemin de réallocation du `vector` s'est écarté du chemin critique. À `-O1`, la machinerie de croissance de `push_back` (la vérification de taille, `operator new`, `memcpy`, `operator delete`) était inlinée en plein milieu de la boucle de collecte. À `-O2` le chemin chaud se réduit à « ajouter le pointeur, incrémenter le curseur, boucler », et toute la séquence de réallocation est déplacée vers un bloc froid en bas, atteint uniquement quand le buffer est réellement plein. Ça ne change pas quand les allocations ont lieu, mais ça garde le cas commun compact. Et `-O3` ? Rien d'intéressant ici. L'autovectorisation ne s'est pas déclenchée (comme on l'a vu, ces boucles ne se vectorisent pas), et il n'y avait aucune boucle qui valait la peine d'être déroulée. Il se retrouve donc juste à côté de `-O2`. ## Performance | Niveau | Frame 0 (à froid) | État stable (~frame 90-270) | Moyenne sur 300 frames | Accélération vs O0 | | ------ | ----------------: | ---------------------------: | ------------------: | --------------: | | **O0** | 1719.47 ms | ~750-790 ms | **778.51 ms** | 1.0x (référence) | | **O1** | 135.35 ms | ~42 ms | **45.13 ms** | 17.3x | | **O2** | 120.82 ms | ~39 ms | **43.90 ms** | 17.7x | | **O3** | 124.03 ms | ~40 ms | **44.22 ms** | 17.6x | Ça correspond à ce que l'assembleur nous disait. La majeure partie du gain vient du saut vers `-O1`, et le reste est minime. Reste que gagner 17x en vitesse sans toucher au code, c'est un bon résultat. ## Mode debug Il existe un flag dans Meson (et dans clang) qui contrôle si le débogage est activé, c'est-à-dire si les symboles de débogage sont inclus. Ça doit sûrement faire une grosse différence, non ? | Run | Frame 0 (à froid) | Moyenne sur 300 frames | Accélération vs O0 | | ------------------- | -------------- | ------------------- | ------------- | | **O3, debug=true** | 124.03 ms | 44.22 ms | 17.6x | | **O3, debug=false** | 135.65 ms | 46.51 ms | 16.7x | Je plaisante. Ce n'est pas le cas. La différence ici n'est que du bruit. Les symboles de débogage se trouvent dans le binaire ; ils ne ralentissent pas le code qui s'exécute. ## Le balayage des flags `-O3` est le levier évident, mais clang a toute une pile d'autres flags d'optimisation. Alors j'ai fait ce que ferait n'importe quel ingé soft : j'ai brute-forcé environ 1025 combinaisons de flags, exécuté chacune d'elles (meilleure sur 3, 300 particules sur 300 frames), et je les ai triées par millisecondes par étape. L'accélération est relative à la référence `-O0` à 47.91 ms/étape. Bien sûr. J'ai demandé à Claude d'écrire un script python pour automatiser tout ça Voici le haut de la liste, plus les lignes `-O3` et `-O0` seules à titre de référence : ``` rank best ms mean ms speedup flags ---------------------------------------------------------------------------------------- 1 1.21 1.23 39.60x -O3 -march=native -ffast-math -funroll-loops -ffp-contract=fast -fno-stack-protector 2 1.21 1.24 39.60x -O3 -march=native -ffast-math -funroll-loops -fno-math-errno -fomit-frame-pointer 3 1.22 1.28 39.27x -O3 -march=native -ffast-math -fno-plt -fno-stack-protector 4 1.22 1.26 39.27x -O3 -march=native -ffast-math -funroll-loops -ffp-contract=fast -fno-math-errno 5 1.22 1.24 39.27x -O3 -march=native -ffast-math -funroll-loops -fno-math-errno -fno-semantic-interposition ... 147 1.26 1.27 38.02x -O3 -march=native -ffast-math ... 945 2.44 2.45 19.64x -O3 (base) ... 1025 47.91 47.99 1.00x -O0 (baseline) ``` La forme de tout ça est ce qu'il y a d'intéressant. `-O3` seul vous amène à environ 19.6x. Tout le haut du tableau se situe autour de 39x. Ce quasi-doublement supplémentaire vient de deux flags, et essentiellement de deux flags seulement : `-march=native` et `-ffast-math`. Remarquez que le rang 147, qui n'est que `-O3 -march=native -ffast-math` sans rien d'autre, atteint déjà 38x. Tout ce qui s'empile par-dessus n'est que de l'erreur d'arrondi. Alors qu'est-ce que tous ces flags font réellement ? à peu près dans l'ordre de leur importance ici. Les deux qui ont porté le résultat : - **`-march=native`** : compiler spécifiquement pour le CPU sur lequel vous compilez, plutôt que pour une base générique. Ça débloque les jeux d'instructions plus récents que votre puce possède réellement (AVX, FMA, et compagnie). Gros gain. - **`-ffast-math`** : permet au compilateur d'assouplir les règles strictes de virgule flottante IEEE. Il peut réassocier les opérations, supposer qu'il n'y a ni NaN ni infinis, et globalement réarranger et fusionner vos calculs de façons que le standard interdit normalement. Ceux qui n'ont presque rien fait sur cette charge de travail : - **`-ffp-contract=fast`** : permet de fusionner une multiplication suivie d'une addition en une seule instruction FMA. Déjà largement couvert par `-ffast-math`. - **`-fno-math-errno`** : empêche les fonctions mathématiques comme `sqrt` et `pow` de positionner `errno`, ce qui permet au compilateur de les traiter comme pures et de les optimiser plus agressivement. - **`-funroll-loops`** : émet plusieurs itérations du corps d'une boucle par passe pour réduire le coût de la boucle. N'a eu que peu d'effet ici, car il n'y avait pas grand-chose qui valait la peine d'être déroulé. - **`-flto`** : optimisation au moment de l'édition de liens, qui permet à l'optimiseur de voir à travers des fichiers source séparés. Négligeable quand la majeure partie du temps cpu se trouve dans une seule boucle d'un seul fichier. - **`-fomit-frame-pointer`** : supprime le frame pointer pour que `rbp` devienne un registre général libre. - **`-fno-plt` et `-fno-semantic-interposition`** : rendent les appels vers des bibliothèques partagées moins coûteux et plus directs. - **`-fno-stack-protector`** : supprime les vérifications du canari de pile. Minime, et vous renoncez à une fonctionnalité de sécurité {{< alert >}} `-ffast-math` et `-fno-stack-protector` ne sont pas gratuits. `-ffast-math` change les résultats réels de vos calculs en virgule flottante, donc si vous vous souciez de la reproductibilité ou de cas limites comme la gestion des NaN, mesurez avant de lui faire confiance. `-fno-stack-protector` supprime une véritable mitigation de sécurité. ASSUREZ-VOUS DE SAVOIR CE QUE VOUS FAITES !!! {{< /alert >}} ## Conclusion Voilà donc la partie facile terminée. Sans réécrire une seule ligne de l'algorithme, les bons flags de compilateur m'ont fait passer de la référence `-O0` à environ 40x plus rapide. La majeure partie vient simplement de `-O3 -march=native -ffast-math`. Mais le compilateur ne peut faire que ce qu'il peut avec le code que je lui ai donné. Il a refusé de vectoriser mes boucles, les données sont éparpillées dans le tas, et la recherche de voisins fait bien plus de travail que nécessaire. C'est là que la partie suivante entre en jeu : rendre le code lui-même rapide.