# Augmenter l'impact de la vectorisation par transformation automatique de programmes

Harris Ludwig Boudi, Philippe Langlois, Guillaume Revy

► To cite this version:

Harris Ludwig Boudi, Philippe Langlois, Guillaume Revy. Augmenter l'impact de la vectorisation par transformation automatique de programmes. ComPAS 2020 - Conférence francophone d'informatique en Parallélisme, Architecture et Système, Jun 2020, Lyon, France. lirmm-04939382

HAL Id: lirmm-04939382

https://hal-lirmm.ccsd.cnrs.fr/lirmm-04939382v1

Submitted on 10 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Augmenter l'impact de la vectorisation par transformation automatique de programmes

Harris Ludwig Boudi, Philippe Langlois et Guillaume Revy

Université de Perpignan Via Domitia, DALI, F-66860, Perpignan, France

Univ. Montpellier, LIRMM, UMR 5506, F-34095, Montpellier, France

CNRS, LIRMM, UMR 5506, F-34095, Montpellier, France

{harris-ludwig.boudi, philippe.langlois, guillaume.revy}@univ-perp.fr

## Résumé

Réduire le format de représentation de certaines données numériques et des instructions d'un programme est une approche de plus en plus répandue pour améliorer ses performances sans trop dégrader sa précision de sortie. Cependant, isoler les parties du programme pour lesquelles cette réduction est pertinente est une tâche compliquée, la combinatoire devenant beaucoup trop importante dès lors que la taille du programme augmente. Cet article présente l'outil VectOptim, qui permet de réduire le format des instructions d'un programme, en favorisant celles qui bénéficieront des unités vectorielles (SIMD). Nous montrons sur l'exemple de la multiplication de matrices qu'un facteur d'accélération de 1.8 peut être atteint automatiquement.

**Mots-clés** : arithmétique virgule flottante, vectorisation, transformation automatique de programme, analyse de performance

## 1. Introduction

Aujourd'hui les applications numériques sont de plus en plus complexes, et la quantité de calcul en arithmétique virgule flottante est alors de plus en plus importante. Écrire des logiciels de calcul scientifique qui réduisent les erreurs numériques tout en tirant le meilleur profit des spécificités d'une certaine micro-architecture est un réel challenge. Par conséquent et pour s'assurer d'une certaine précision sur le résultat, ces logiciels sont très souvent écrits en utilisant la plus grande précision de calcul directement disponible en matériel, ce indépendamment de la précision réellement souhaitée, choix effectué au détriment de la performance.

Pour améliorer les performances de ces logiciels sans trop dégrader la précision fournie, une approche de plus en plus répandue est d'utiliser la précision mixte. Un programme en précision mixte est un programme qui combine des données dans différents formats flottants, typiquement les formats binary32 et binary64, voire binary16 [1]. L'intérêt d'une telle approche est double. (1) Elle permet de limiter le coût dû au chargement de certaines données de la mémoire. (2) Elle permet également d'exploiter au mieux les spécificités des micro-architectures modernes, notamment le support aux instructions SIMD. Utilisée "à la main", dans le cas de la résolution de systèmes linéaires denses par raffinement itératif, cette approche a par exemple permis d'observer une accélération d'un facteur 1.8 entre un programme en précision mixte binary32/64 et un en binary64 [3].
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Opération scalaire vs. opération vectorielle, en AVX2.</figcaption></figure>

Mais la principale difficulté pour automatiser ce type d'approche vient de la combinatoire des possibilités (e.g. Pⁿ possibilités, avec P formats et n données). Récemment de nombreux outils automatiques qui permettent d'adapter le format des données d'un programme pour satisfaire des contraintes de précision ou de performance ont été développés. Nous les présentons en Section 2. Dans cet article, l'approche proposée pour diminuer la combinatoire est de ne s'intéresser qu'aux *nids de vectorisation*, i.e. aux instructions qui bénéficieront des unités SIMD. Implantée dans l'infrastructure de compilation LLVM, elle permet d'augmenter l'impact de la vectorisation et d'atteindre, par exemple, un facteur 1.8 sur la multiplication de matrices.

L'article est organisé comme suit : La Section 2 rappelle quelques éléments de l'état de l'art. La Section 3 présente l'outil d'optimisation VectOptim. Enfin, la Section 4 montre quelques premiers résultats expérimentaux, avant de conclure en Section 5.

## 2. État de l'art

**Arithmétique à virgule flottante IEEE 754 et SIMD.** L'arithmétique à virgule flottante est le système de représentation le plus largement utilisé pour manipuler les nombres réels en informatique. Elle est décrite dans le standard IEEE 754-2008 [1]. Un nombre flottant binaire x est donné par un triplet (s, e, m) avec $x = (-1)^s \cdot 2^e \cdot m$, où s ∈ {0, 1} est le signe, e ∈ [e_min, e_max] est un exposant entier relatif, et m est la mantisse binaire codée sur p bits. L'intervalle [e_min, e_max] et la taille p de la mantisse dépendent du format de représentation utilisé. Le standard définit différents formats, dont les formats binary32 (simple précision) et binary64 (double précision), que nous considérons dans cet article. Sans entrer dans le détail, notons que les données dans ces deux formats sont respectivement encodés sur 32 et 64 bits.<fnref n="1" />

Ici, nous utilisons l'existence de ces deux représentations pour mieux tirer profit des instructions SIMD, et notamment augmenter leur impact sur les performances des applications en réduisant le format de certaines instructions. En effet le SIMD (Single Instruction Multiple Data) est un paradigme qui exploite le parallélisme des applications en exécutant de façon simultanée une instruction sur plusieurs données. Ces données sont organisées en vecteurs, dont la taille varie en fonction du jeu d'instructions utilisés, e.g. SSE, AVX, AVX2, ou AVX512 pour Intel. Par exemple, les micro-architectures supportant un jeu d'instructions AVX2 ont des registres vectoriels de 256 bits, qui permettent de manipuler 4 binary64 ou 8 binary32 en une seule instruction. Donc intuitivement, réduire le format des instructions binary64 vers le format binary32, i.e. réduire d'un facteur 2 la taille des données manipulées, permet d'exécuter le double d'opérations élémentaires à chaque instruction, comme illustré sur la Figure 1, et ainsi d'améliorer d'un facteur 2 les performances de certaines parties du programme.

**Principaux outils existants.** La réduction automatique de formats flottants est un domaine de recherche très actuels. Plusieurs approches ou outils ont été récemment publiés. Craft [9], intégré dans le projet FloatSmith,<fnref n="2" /> permet, par instrumentation binaire, d'isoler des parties d'un

<footnote n="1">1. Plus de détails sur les formats de données et leur encodage sont dans le standard IEEE 754-2008 [1].</footnote>

<footnote n="2">2. Informations disponibles à l'adresse https://github.com/crafthpc/floatsmith.</footnote>
programme dont la précision peut être réduite. Il explore différentes configurations par une recherche en profondeur d'abord. Precimonious [14] et Blame [13] sont deux outils qui utilisent un algorithme de delta-debugging [16] pour déterminer quelles données peuvent être transformées en précision réduite. Pour identifier quelle configuration est la plus pertinente, ces outils comparent simplement chaque résultat en précision mixte au résultat en précision double. Promise [6] utilise aussi un algorithme de delta-debugging mais repose sur l'arithmétique stochastique pour évaluer l'impact des transformations (via l'outil CADNA [7]). Citons AMP [12] qui, contrairement aux précédents, augmente la précision de certaines portions du programme pour améliorer la précision du résultat sans impacter les performances. Une autre approche est implantée dans STOKE [15], qui applique des transformations aléatoires, déterminées par des méthodes stochastiques. Plus récemment, AMPT-GA a été proposé [8] : il utilise des algorithmes génétiques pour minimiser une fonction de l'erreur et des performances sur GPU. Ces outils fonctionnent de manière dynamique, en exécutant le programme sur un ou plusieurs jeux de données. Une autre approche consiste à travailler au niveau du source : FP-Tuner [5] utilise de l'analyse formelle via les séries de Taylor symboliques, alors que dans [11], des contraintes de précision sur les données sont déterminées par analyse statique, pouvant en-suite être vérifiées par un solveur SMT. DAISY utilise également des solveurs SMT, combinés avec de l'arithmétique affine et par intervalle [4]. Le principale inconvénient de ces techniques est qu'elles ne permettent pas de réellement réduire à la combinatoire des possibilités, et elles deviennent coûteuses, dès que la taille du programme augmente.

### 3. VectOptim

Nous présentons maintenant le projet VectOptim. Il intègre un flot de transformations automatiques de programmes, qui réduit le format des instructions apparaissant dans les *nids de vectorisation*. Implanté dans LLVM [10] comme une *passe* de transformation, il fonctionne en trois étapes. Après une présentation des capacités offertes par l'outil LLVM, cette section présente ces trois étapes et leurs interactions.

#### 3.1. Infrastructure LLVM

LLVM (Low Level Virtual Machine) propose une infrastructure de compilation particulièrement adaptée à l'analyse, la transformation et l'optimisation de programmes. Son fonctionnement suit celui d'un compilateur classique [2], en trois étapes : le front-end (i.e. clang pour les langages C et C++), le middle-end (en réalité, LLVM), et le back-end. LLVM repose sur une représentation intermédiaire en forme SSA (Static Single Assignment), i.e. chaque variable (appelée *registre typé*) est affectée une seule fois. Obtenue en sortie du front-end, cette représentation utilise un jeu d'instructions fortement typé et indépendant du langage de programmation et de l'architecture cible, donc particulièrement adaptée à l'optimisation de programmes.

Comme la plupart des compilateurs modernes, LLVM intègre un certain nombre de phases d'analyse et de transformation, incluant celles d'optimisation classique, comme la vectorisation de boucles. Dans LLVM, ces phases sont implantées sous formes de *passes*. Mais ce qui différencie LLVM de ses principaux concurrents, notamment GCC, est le mécanisme qu'il propose pour le développement de nouvelles passes,<fnref n="3" /> mécanisme que nous utilisons de ce travail.

#### 3.2. Identification des nids de vectorisation

La première étape de VectOptim consiste à identifier les *nids de vectorisation*. En effet, la plupart des compilateurs actuels intègrent une phase d'auto-vectorisation permettant d'adapter au-

<footnote n="3">3. Plus de détails à l'adresse http://llvm.org/docs/WritingAnLLVMPass.html.</footnote>
tomatiquement un programme scalaire à un processeur SIMD. Les nids de vectorisation sont alors des ensembles d'instructions scalaires qui, une fois cette phase de vectorisation effectuée, seront transformées en instructions vectorielles.

Ces nids de vectorisation sont facilement identifiables sur une représentation intermédiaire vectorisée : il suffit d'identifier toutes les instructions vectorielles. Mais il est plus efficace d'appliquer nos transformations sur le code scalaire. <fnref n="4" /> Nous commençons par vectoriser le programme initial, et nous utilisons ensuite le mécanisme de *méta-données* de LLVM pour remonter des instructions vectorielles aux nids de vectorisation correspondants dans le programme scalaire. Ce mécanisme permet de marquer les instructions de la représentation intermédiaire avec un tag. Considérons le code C ci-dessous (à gauche), et un aperçu de la représentation intermédiaire *marquée* (à droite).

```c
for (int i = 0; i < n; i++)
    c[i] = a[i] + b[i];
```

```c
for body:                      ; preds = %for.cond
    ; construction des operandes %6 (a[i]) et %9 (b[i])
    %fadd.0 = fadd double %6, %9, !vectoptim.orig.name !50
    ; construction de l'opérande %arrayidx6 (c[i])
    store double %fadd.0, double* %arrayidx6, align 8
    br label %for inc
// ...
!50 = !!"fadd.0"
```

Dans cet exemple, l'instruction d'addition est renommée de manière unique (ici, `fadd.0`), puis marquée avec un tag, dont le nom est `vectoptim.orig.name` (en bleu) et dont la valeur (`!50`, en rouge) fait référence au nom de l'instruction dans la représentation intermédiaire scalaire (en vert). Une fois la vectorisation effectuée (e.g. pour un jeu d'instructions AVX2), cette instruction sera transformée en une instruction vectorielle avec le même tag, comme illustré ci-dessous.

```c
vector.body:                      ; preds = %vector.body, %vector.ph
    %34 = fadd <4 x double> %22, %wide.load37, !vectoptim.orig.name !43
// ...
!43 = !!"fadd.0"
```

Nous pouvons alors déduire que l'instruction d'addition `%34` est issue de l'addition `fadd.0`, qui fait partie d'un nid de vectorisation : elle peut donc être transformée.

Notons tout de même que nous ne marquons bien évidemment pas toutes les instructions scalaires, mais seulement celles qui nous intéressent, à savoir certaines instructions flottantes (opérations arithmétiques et appels aux fonctions mathématiques, notamment) au format binary64.

### 3.3. Réduction du format des instructions

La deuxième étape consiste à réduire le format des instructions qui apparaissent dans les nids de vectorisation précédemment identifiés. Par réduire le format, nous entendons les remplacer par des instructions dans un format plus court. En effet, nous avons vu en Section 3.1 que LLVM utilise une représentation intermédiaire fortement typée, ce qui signifie que chaque instruction possède son propre format. Mais la documentation de LLVM indique que modifier le type d'un objet (variable, instruction, ...) peut être une "opération extrêmement dangereuse, qui peut produire une représentation intermédiaire complètement invalide", et qu'il est "fortement recommandé de re créer un nouvel objet avec le type correct". <fnref n="5" /> Ainsi, au lieu de modifier le type d'une instruction de binary64 vers binary32, nous dupliquons cette instruction en son équivalent en format réduit, que nous insérons ensuite dans la représentation intermédiaire, comme illustré sur la Figure 2. Notons que nous insérons également les opérations de conversions, binary64 vers binary32 (→b32) et binary64 vers binary32 (→b32), pour garantir une cohérence entre les instructions et leurs opérandes, en prenant soin de supprimer les conversions opposées qui se suivent, comme illustré en bleu sur la Figure 2(b).

<footnote n="4">4. Du moins, sur une représentation intermédiaire correspondant au code scalaire, sans aucune optimisation.</footnote>

<footnote n="5">5. Voir la méthode `mutateType` à l'adresse https://llvm.org/doxygen/classllvm_1_1Value.html.</footnote>
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Méthode pour la réduction du format des instructions flottantes.</figcaption></figure>

## 3.4. Amélioration de la vectorisation par duplication de données

La troisième étape consiste à dupliquer certaines données pour favoriser la vectorisation. En effet, nous pouvons observer que les nids de vectorisation sont principalement des ensembles d'instructions situées dans les corps de boucles. Les transformations introduites précédemment impliquent une conversion $\rightarrow_{b32}$ pour chaque opérande et une conversion $\rightarrow_{b64}$ pour convertir le résultat de la nouvelle instruction, comme illustré sur le code C ci-dessous (à gauche).

```c
for (int i = 0; i < n; i++)
    a32 = (float)a[i]; // conversion $\rightarrow_{b32}$
b32 = (float)b[i]; // conversion $\rightarrow_{b32}$
c32 = a32 + b32;
c[i] = (double)c32[i]; // conversion $\rightarrow_{b64}$
```

```c
COPY(a, a32); // conversion $\rightarrow_{b32}$ pour a
COPY(b, b32); // conversion $\rightarrow_{b32}$ pour b
for (int i = 0; i < n; i++)
    c32[i] = a32[i] + b32[i];
COPY(c32, c); // conversion $\rightarrow_{b64}$ pour c
```

Ces conversions sont évaluées à chaque itération de boucle, ce qui peut être pénalisant, à deux niveaux. (1) Ce sont des opérations coûteuses, qui peuvent être exécutées plusieurs fois de manière redondante, et donc impacter négativement les performances. (2) Elles peuvent également empêcher de bénéficier pleinement des unités vectorielles. Plus précisément, dans la représentation intermédiaire de LLVM, les conversions $\rightarrow_{b32}$ sont des opérations binary64. Or, pour chaque instruction de la boucle, le vectoriseur détermine un facteur de vectorisation, indiquant la taille des vecteurs qu'elle peut manipuler. Par exemple, en AVX2, dans le cas d'une conversion $\rightarrow_{b32}$, le facteur de vectorisation est de 4 (i.e. $256 = 4 \times 64$ bits). Ce facteur étant au final commun pour toutes les instructions d'une même boucle, il est également de 4 pour les instructions binary32, alors qu'il pourrait être de 8 (i.e. $256 = 8 \times 32$ bits). Ceci empêche de bénéficier des instructions "$8 \times$ binary32". Pour remédier à cela, nous déplaçons ces conversions en dehors des boucles, en s'autorisant à dupliquer certaines données. Par exemple, sur le programme C ci-dessus (à droite), nous déportons les conversions $\rightarrow_{b32}$ avant le corps de boucle, et les conversions $\rightarrow_{b64}$ après pour mettre à jour les éléments $c[i]$. Ainsi les conversions sont effectuées une fois pour toutes, et n'impactent par la vectorisation, d'autant plus qu'elles font apparaître elles-mêmes des boucles vectorisables.

## 4. Premiers résultats expérimentaux

Nous avons implanté cette démarche dans l'infrastructure LLVM 8.0, et nous présentons ici son utilisation sur l'exemple de la multiplication de matrices $C = A \times B$, avec A, B et C de taille $n \times n$ à coefficients dans le format binary64. Dans cette expérience, nous considérons 5 variantes du programme. Les 3 premières reposent sur l'implantation "naïve", avec trois boucles imbriquées : la version originale en binary64 ($V_1$), la version optimisée par notre outil ($V_2$), et la version en binary32 ($V_3$). Les 2 dernières utilisent la méthode par bloc, avec un découpage
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Speedup et erreur relative par rapport à la version originale en binary64 ($V_1$).</figcaption></figure>

en 4 blocs : une version optimisée ($V_4$), et une version où seule la multiplication des blocs est optimisée, l'addition étant effectuée en binary64 ($V_5$). Chaque programme est compilé en niveau 3 d'optimisation, avec support AVX2 activé, et exécuté sur un ensemble de 1000 couples de matrices (A, B) aléatoires. En sortie, nous mesurons le speedup (ratio entre le nombre de cycles minimum mesuré avec RDTSC en binary64 et celui pour le programme transformé), et l'erreur relative maximum sur les éléments de C (par rapport à la matrice résultat calculée en binary64). La Figure 3 montre ces mesures pour des matrices de petite taille variant de 16 à 256. Nous observons que, malgré la duplication des données, tous les programmes transformés sont plus rapides que la version originale en binary64 ($V_1$), excepté la multiplication par bloc pour $n \le 64$, voire parfois même plus rapide que la version binary32 ($V_3$). Notamment, pour la version optimisée ($V_2$), un speedup de $\approx 1.8$ peut être mesuré, ce qui semble pertinent dans notre contexte. Concernant la multiplication par bloc, la version optimisée où seule la multiplication est en binary32 ($V_5$) est certes plus rapide que la version binary64 ($V_1$), mais surtout, de manière contre-intuitive, plus rapide que la version par bloc complètement optimisée ($V_4$). Ceci renforce l'idée qu'il n'est pas forcément nécessaire de transformer la totalité du programme pour améliorer ses performances. D'un point de vue précision, la version optimisée ($V_2$) est, comme attendu, équivalente à la version binary32 ($V_3$). Quant aux deux versions par bloc, elles restent sensiblement équivalentes sur l'ensemble de nos tests. Dans ce sens, la version $V_5$ est plus rapide et aussi précise, donc plus pertinente que la version $V_4$. Cependant, décider quelles portions du programme doivent ou non être transformées est fait ici "à la main", ce qui peut être, sur un programme plus conséquent, fastidieux. Ceci justifie encore, si il était besoin, l'intérêt d'outils automatiques comme ceux étudiés ici.

## 5. Conclusion

Cet article présente VectOptim, outil qui permet de réduire le format des instructions flottantes d'un programme pouvant bénéficier des unités SIMD. Sur l'exemple de la multiplication de matrices, nous observons un gain intéressant d'un facteur 1.8 par rapport à la version initiale. Des travaux futurs s'envisagent à trois niveaux. (1) Nous devrons évaluer cet outil sur des programmes plus conséquents pour observer le gain qu'il permet réellement. (2) Notons que ce travail est une première étape vers un outil d'adaptation automatique des formats de données, et qu'il sera intéressant de l'intégrer dans un flot complet de transformations qui intègre la satisfaction des contraintes de précision du programme ainsi transformé. (3) Finalement, nous nous sommes intéressé à la transformation binary64 vers binary32. Il serait certainement pertinent d'étudier également la possibilité de transformer certaines instructions vers le binary16, et ainsi atteindre un gain encore plus conséquent.
Bibliographie

1. IEEE standard for floating-point arithmetic. - IEEE Std. 754-2008, pages 1-58, 2008.

2. Aho (A. V.), Sethi (R.) et Ullman (J. D.). – *Compilers: principles, techniques and tools*. – Addison Wesley, 1986.

3. Baboulin (M.), Buttari (A.), Dongarra (J.), Kurzak (J.), Langou (J.), Langou (J.), Luszczek (P.) et Tomov (S.). – Accelerating scientific computations with mixed precision algorithms. *Computer Physics Communications*, vol. 180, n12, 2009, pp. 2526–2533.

4. Becker (H.), Darulova (E.) et Myreen (M. O.). – A verified certificate checker for floating-point error bounds. 2017.

5. Chiang (W.-F.), Baranowski (M.), Briggs (I.), Solovyev (A.), Gopalakrishnan (G.) et Raka-marić (Z.). – Rigorous floating-point mixed-precision tuning. – In *POPL 2017*, pp. 300–315, 2017.

6. Graillat. (S.), Jézéquel (F.), Picot (R.), Févotte (F.) et Lathuilière (B.). – Auto-tuning for floating-point precision with Discrete Stochastic Arithmetic. – 2016.

7. Jézéquel (F.) et Chesneaux (J.-M.). – CADNA : a library for estimating round-off error propagation. *Computer Physics Communications*, 2008, pp. 933–955.

8. Kotipalli (P.), Singh (R.), Wood (P.), Laguna (I.) et Bagchi (S.). – Ampt-ga : automatic mixed precision floating point tuning for gpu applications. – 06 2019.

9. Lam (M. O.), Hollingsworth (J. K.), de Supinski (B. R.) et Legendre (M. P.). – Automatically adapting programs for mixed-precision floating-point computation. – In *ICS'13*, 2013.

10. Lattner (C.) et Adve (V.). – LLVM : a compilation framework for lifelong program analysis and transformation. – In *International Symposium on Code Generation and Optimization*, 2004. CGO 2004., pp. 75–86, 2004.

11. Martel (M.). – Floating-Point Format Inference in Mixed-Precision. – In *9th NASA Formal Methods Symposium.*, 2017.

12. Nathan (R.), Naeimi (H.), Sorin (D. J.) et Sun (X.). – Profile-driven automated mixed precision. *CoRR*, vol. abs/1606.00251, 2016.

13. Rubio-González (C.), Nguyen (C.), Mehne (B.), Sen (K.), Demmel (J.), Kahan (W.), Iancu (C.), Lavrijsen (W.), Bailey (D. H.) et Hough (D.). – Floating-point precision tuning using Blame analysis. – In *ICSE 2016*, pp. 1074–1085, 2016.

14. Rubio-González (C.), Nguyen (C.), Nguyen (H. D.), Demmel (J.), Kahan (W.), Sen (K.), Bailey (D. H.), Iancu (C.) et Hough (D.). – Precimonious : tuning assistant for floating-point precision. – In *SC'13*, pp. 1–12, 2013.

15. Schkufza (E.), Sharma (R.) et Aiken (A.). – Stochastic Optimization of Floating-point Programs with Tunable Precision. – In *PLDI'14*, pp. 53–64, 2014.

16. Zeller (A.) et Hildebrandt (R.). – Simplifying and Isolating Failure-Inducing Input. *IEEE Trans. Softw. Eng.*, 2002, pp. 183–200.