Optimiser une simulation SPH, partie 1 : la vitesse offerte par le compilateur

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.

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 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 :

lea     rax, [rbp-0x40]
mov     rsi, rcx
mov     rdi, rax
call    glm::operator-<double, 0>
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 :

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() :

call    std::vector<sph::Particle>::begin
mov     qword [rbp-0x1e0 {__for_begin}], rax
call    std::vector<sph::Particle>::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 :

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 #

Ç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 ?

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é

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.