Les choses sérieuses commencent. Maintenant que Phi n’a plus de secrets pour nous, il est temps de solliciter ses ressources avec, ce mois-ci, un long focus sur la vectorisation et son optimisation.

Dans notre numéro précédent, nous avons commencé notre prise en mains de Phi par un tour complet du propriétaire, en nous attardant sur les ressources matérielles propres à l’accélérateur et sur ce qu’elles impliquent du point de vue de la programmation parallèle. Ce mois-ci, nous allons examiner plusieurs voies d’optimisation liées plus particulièrement à la vectorisation. Considérons-les comme autant d’outils élémentaires – et de réflexes à adopter – quel que soit le type d’application sur laquelle vous travailliez.

Vous le savez probablement déjà, les registres et les opérations SIMD sur 512 bits sont un des apports majeurs de Phi pour le développeur x86. Si l’on veut en tirer le meilleur profit, la première étape consiste à prendre l’habitude de coder en utilisant des notations de type tableau, comme on le fait en Fortran 90 ou avec Intel Cilk Plus. Cette notation garantit que le compilateur vectorisera les instructions en utilisant le jeu d’instruction SIMD partout où ce sera possible.

Petit rappel, la syntaxe de base pour le référencement d’un tableau ou d’un segment de tableau est la suivante :

[<limite inférieure> : <longueur> : <alignement>];

Dans cette définition, <limite inférieure> est le premier élément du tableau sur lequel on opère, <longueur> est le nombre d’éléments du tableau sur lesquels on opère et <alignement> est la distance entre chaque élément du tableau (soit l’équivalent de “stride” en anglais) – valeur généralement égale à 1. Si l’on opère sur l’ensemble du tableau, la spécification des valeurs peut être omise.  Ainsi, les deux lignes suivantes sont équivalentes :

 A[:] = c * B[:];  A[0:n:1] = c * B[0:n:1]; 

et correspondent à la notation classique :

 for (i=0; i<n; i++)    A[i] = c * B[i]; 

Elément par élément

Passons maintenant aux fonctions élémentales. Les extensions Intel pour C et C++ nous permettent de définir nos propres fonctions élémentales, comme on le fait naturellement en Fortran. Pour mémoire, une fonction est dite élémentale si elle peut être invoquée de façon parallèle sur des éléments scalaires ou appartenant à des tableaux. Lorsque l’on déclare une fonction élémentale, le compilateur génère deux versions de la fonction : une première version scalaire normale et une seconde dans laquelle les données sont traitées en parallèle. C’est cette seconde version qui est appelée depuis les boucles for ou lorsque les données à traitées sont vectorielles. Pour définir une fonction élémentale, il faut spécifier l’attribut vector dans sa déclaration, comme dans l’exemple suivant :

__attribute__(vector (optional clauses))  void MyVecMult(double *a, double *b, double *c) {    c[0] = a[0] * b[0];   return; }

Une fois déclarée, notre fonction pourra être appelée comme ceci :

 for (i=0; i<n; i++)   MyvecMult(a[i], b[i], c[i]); 

ou encore comme cela :

 MyvecMult(a[:], b[:], c[:]); 

Notez que des clauses supplémentaires peuvent être ajoutées à la déclaration pour spécifier des longueurs de vecteurs prédéfinies et d’autres “détails” de nature à augmenter la lisibilité du code ou à aider le compilateur à prendre les bonnes décisions. Notez également – et peut être plus encore – qu’il existe un certain nombre de restrictions à ce qu’une fonction élémentale peut accomplir : pas d’instructions switch, pas de goto… Pour connaître la liste complète de ces interdits et des raisons qui les motivent, référez-vous à un excellent article (en anglais) publié sur le site d’Intel.

Domestiquer le compilateur

A ce stade, vous pouvez vous demander si cette notation “tableau” est absolument indispensable pour bénéficier de la vectorisation automatique. En théorie, la réponse est non. Mais en pratique, faire aveuglément confiance aux facultés d’auto-vectorisation du compilateur n’est vraiment recommandé. La bonne pratique, comme on le laissait entendre plus haut, est d’aider le compilateur aussi souvent que c’est possible – c’est-à-dire de prendre le plus grand nombre de décisions vous-même et de les lui indiquer en utilisant les pragmas ou les directives appropriées.

Avant, donc, de passer à l’optimisation de votre code, profilez l’application (si vous utilisez Xeon Phi, vous utilisez sans doute aussi VTune Amplifier XE), activez les rapports de vectorisation à la compilation (via, par exemple, -vecreport3) et vérifiez que les hotspots qui vous sont indiqués sont bien vectorisés. S’ils ne le sont pas, c’est que le compilateur doute : donc, une réécriture en notation tableau ou l’ajout d’un certain nombre de directives et de pragmas s’impose.

Par exemple, si un hotspot avec une boucle for n’est pas vectorisé et si vous préférez l’option déclarative, l’ajout de #pragma ivdep demandera explicitement au compilateur d’ignorer les dépendances de pointeurs potentielles ou supposées (ivdep signifiant littéralement ignore vector dependence). Bien sûr, il faudra vous assurer préalablement que les pointeurs en question déréférencent toujours des zones mémoire indépendantes. Une autre approche consiste à utiliser #pragma simd. Cette directive est moins dirigiste : elle demande au compilateur d’ignorer tous les conflits et de produire du code mobilisant les instructions SIMD partout où il le peut. Là encore, comme c’est vous, le développeur, qui demandez au compilateur d’ignorer tout conflit qu’il ne pourrait pas désambiguer, c’est à vous de veiller à ce que ne subsiste aucune dépendance possible dans le code.

Optimiser la vectorisation

Une fois prises ces précautions de base, on peut passer à l’étape suivante : vérifier que la vectorisation produite est la plus efficace possible. La règle, vous vous en doutez, est que plus le compilateur dispose d’informations sur le code à traiter, plus il pourra travailler en finesse. Pour commencer, regardez si les accès mémoire (en lecture et en écriture) sont bien regroupés et vérifiez que le stockage temporaire des données nécessite le moins possible d’indexations (à moins, bien sûr, que ces indexations soient indispensables comme dans le cas de nombreuses opérations matricielles).

L’alignement des vecteurs est lui aussi extrêmement important pour la performance globale de vos applications. Quand vous déclarez des tableaux, veillez à les aligner sur des adresses 64 bits, en utilisant une déclaration du type __declspec(align(64)) float A[1000]; ou en faisant appel, en C, à l’instruction _aligned_malloc(). Lorsque vous échangez des pointeurs entre fonctions ou routines, pensez à utiliser la #pragma vector aligned immédiatement avant la boucle, pour indiquer au compilateur que tous les pointeurs sont effectivement alignés. S’ils le sont pour toutes les boucles d’une routine, vous pouvez aussi déclarer __assume_aligned() pour le pointeur, ce qui évite la déclaration précédente pour chaque boucle. Enfin, lorsque deux boucles imbriquées se présentent et qu’elles sont courtes, essayez de les regrouper en une seule. Pour les raisons évoquées plus haut, ce petit travail manuel est préférable à l’utilisation de pragmas telles que unroll_and_jam et nounroll_and_jam.

N’oublions pas la mémoire

Dernier grand domaine général d’optimisation, la gestion de la mémoire. Clairement, les meilleures conditions d’exécution du code sont celles où les accès aux données se font sur des plages d’adresses séquentielles. Pour les développeurs Phi les plus expérimentés, restructurer les données dans le but d’obtenir cette séquentialité est une activité fréquente. Généralement, elle consiste à passer d’un tableau de structures à une structure de tableaux.

Sachant par ailleurs que le prefetching fait beaucoup pour la réduction du temps d’exécution, il faut le soigner. La plupart des compilateurs assurent automatiquement un prefetching sur les boucles qu’ils vectorisent. Mais si le pattern d’accès mémoire de telle ou telle boucle n’est pas clair à la compilation, il peut être intéressant de spécifier explicitement ce prefetch – par pragma, par directive ou par commande intrinsèque. A ce sujet, n’oubliez pas que si vous insérez une instructions de prefetch dans une boucle for, cette instruction doit viser la prochaine itération, pas l’itération courante, sinon elle ne donnera aucun résultat effectif.

Pour affiner encore les choses, essayez autant que faire se peut de réutiliser les données pendant qu’elles sont encore dans les registres ou dans le cache du processeur. Il est très souvent possible de partitionner les boucles de façon à ce que le cache n’ait pas besoin d’éjecter les données dont on a un besoin immédiat, et de grouper ces mêmes données en blocs de la taille du cache local concerné (512 Ko par exemple). Dans l’exemple ci-dessous, en notation tableau, la seconde implémentation est ainsi sensiblement optimisée :

 A[0:N] = B[0:N]+C[0:N]; D[0:N] = E[0:N]+A[0:N];  #define VLEN 4 for(i=0;i<N;i+=VLEN){   A[i:VLEN] = B[i:VLEN]+C[i:VLEN];   D[i:VLEN] = E[i:VLEN]+A[i:VLEN]; }

Bons développements !

Trois cas typiques d’optimisation

1 – Gérer la désambiguation mémoire à l’intérieur des boucles vectorisées

Prenons un exemple de vectorisation sur une boucle simple :

 void vectorize(float *a, float *b, float *c, float *d, int n) {   int i;   for (i=0; i<n; i++) {     a[i] = c[i] * d[i];     b[i] = a[i] + c[i] - d[i];   } }

Dans cette implémentation, le compilateur n’a aucune idée de ce vers quoi les quatre pointeurs pointent. Intuitivement, le développeur sait qu’ils pointent vers des adresses complètement différentes. Le compilateur, en revanche, ne peut pas le deviner. Par conséquent, il présuppose que ces pointeurs sont mal aliasés – que par exemple c[1] et a[0] peuvent être à la même adresse, et que de ce fait la boucle ne peut pas être vectorisée.

Lorsque que le nombre de pointeurs est très réduit (moins d’une dizaine), le compilateur peut générer un code de vérification à l’exécution et deux versions de la boucle (optimisée et non-optimisée). Evidemment, cette approche se révèle coûteuse en temps et en taille de code, de sorte qu’il faut mieux faire l’effort d’indiquer au compilateur que les pointeurs sont indépendants.

Pour cela, utilisez le mot-clé restrict pointer de C99. Et si vous ne compilez pas en “C99 standard“, vous pouvez toujours utiliser les flags -restrict (Linux) et -Qrestrict (Windows) pour que le compilateur accepte le mot-clé restrict. Ainsi, pour spécifier que a et b ne sont aliasés avec rien d’autres, il suffit d’écrire :

 void vectorize(float *restrict a, float *restrict b, float *c, float *d, int n) {   int i;   for (i=0; i<n; i++) {     a[i] = c[i] * d[i];     b[i] = a[i] + c[i] - d[i];   } }

Vous pouvez également utiliser la pragma ivdep. La sémantique est légèrement différente par rapport à restrict mais permet néanmoins au compilateur d’éliminer les assomptions de dépendances et donc de considérer que la vectorisation est possible :

 void vectorize(float *a, float *b, float *c, float *d, int n) {   int i;   #pragma ivdep   for (i=0; i<n; i++) {     a[i] = c[i] * d[i];     b[i] = a[i] + c[i] - d[i];   } }

2 – Vectoriser des boucles avec accès indirect

Prenons l’exemple d’une boucle implémentant des accès mémoire (lecture / écriture) indirects :

 for (i = kstart; i < kend; ++i) {   istart = iend;   iend   = mp_dofStart[i+1];   float w = xd[i];    for (j = istart; j < iend; ++j) {     index  = SCS[j];     xd[index] -= lower[j]*w;   } }

Dans cette implémentation, le pré-requis majeur pour la vectorisation est que les valeurs xd soient distinctes. Si elles ne le sont pas, les dépendances empêcheront la vectorisation de la boucle. Si elles le sont, vous pouvez utiliser les pragmas ivdep/simd (avant la boucle j), qui donneront une vectorisation plus efficace à la compilation.

3 – Vectoriser des boucles avec variable d’induction monotonique

Les compilateurs peuvent vectoriser certaines boucles utilisant une variable d’induction monotonique (c’est-à-dire une variable qui n’est mise à jour que sous certaines conditions d’exécution de la boucle). En voici un exemple :

 int index_0 = 0; for(int k0=0; k0 < count0; k0++) {   TYPE X1 = *(Pos0 + k0);           TYPE Y1 = *(Pos0 + k0 + count0);   TYPE Z1 = *(Pos0 + k0 + 2*count0); #pragma loop_count min(220) avg (300) max (380) #pragma ivdep   for(int k1=0; k1 < count1; k1+=1) {     TYPE X0 = *(Pos1 + k1);     TYPE Y0 = *(Pos1 + k1 + count1);      TYPE Z0 = *(Pos1 + k1 + 2*count1);     TYPE diff_X = (X0 - X1);     TYPE diff_Y = (Y0 - Y1);     TYPE diff_Z = (Z0 - Z1);     TYPE norm_2 = (diff_X*diff_X) + (diff_Y*diff_Y) + (diff_Z*diff_Z);           if ((norm_2 >= rmin_2) && (norm_2 <= rmax_2))       Packed[index_0++] = norm_2;   } }

Pour ce type de boucle, on peut améliorer la génération de code via l'option -opt-assume-safe-padding. Quand elle est spécifiée, le compilateur part du principe que les variables et la mémoire allouée dynamiquement sont paddées au-delà des limites de l'objet. Il permet donc au code d'accéder à 64 octets de plus que ce que le programme indique partout où l'objet est présent. Attention toutefois, pour que tout se passe bien, il est nécessaire d'augmenter la taille des objets statiques et dynamiques.

Quel est le bénéfice concret de l'opération ? Prenons l'exemple d'une boucle similaire à celle qui illustrait le bénéfice des pointeurs restrict :

 void foo(float* restrict a, float* restrict b, float* restrict c) {   int i;   int j = 0;   for(i = 0; i < N; i++) {     if (b[i])  {       a[j++] = c[i];     }   } }

Sans l'option -opt-assume-safe-padding (c'est-à-dire en configuration par défaut), le compilateur se montre prudent (pour éviter les erreurs mémoire) et génère le code vectorisé suivant :

 ..B1.6:   vloadunpackld (%rsi,%rax,4), %zmm2   vprefetch1 512(%rsi,%rax,4)   vloadunpackld (%rdx,%rax,4), %zmm3   vprefetch0 256(%rsi,%rax,4)   vloadunpackhd 64(%rsi,%rax,4), %zmm2   vprefetch1 512(%rdx,%rax,4)   vloadunpackhd 64(%rdx,%rax,4), %zmm3   vprefetch0 256(%rdx,%rax,4)   vcmpneqps %zmm1, %zmm2, %k1   movl      $65535, %r10d   vpackstorelps %zmm3, -64(%rsp){%k1}   lea       (%rdi,%r8,4), %r11   vmovaps   -64(%rsp), %zmm4   kmov      %k1, %r9d   popcnt    %r9d, %ecx   addq      $16, %rax   movl      %ecx, %ecx   shll      %cl, %r10d   addq      %rcx, %r8   xorl      $-1, %r10d   kmov      %r10d, %k2 ..L7:   vscatterdps %zmm4, (%r11,%zmm0,4){%k2}   jkzd      ..L6, %k2   vscatterdps %zmm4, (%r11,%zmm0,4){%k2}   jknzd     ..L7, %k2 ..L6:   cmpq      $992, %rax   jb        ..B1.6 

Avec l'option -opt-assume-safe-padding (et en partant du principe que vous avez ajouté 64 octets de padding lors de l'allocation du tableau), le code généré par le compilateur devient nettement plus performant :

 ..B1.6:   vloadunpackld (%rsi,%rax,4), %zmm1   vprefetch1 512(%rsi,%rax,4)   vloadunpackld (%rdx,%rax,4), %zmm2   vprefetch0 256(%rsi,%rax,4)   vloadunpackhd 64(%rsi,%rax,4), %zmm1   vprefetch1 512(%rdx,%rax,4)   vloadunpackhd 64(%rdx,%rax,4), %zmm2   vprefetch0 256(%rdx,%rax,4)   vcmpneqps %zmm0, %zmm1, %k1   movl      $65535, %r10d   vpackstorelps %zmm2, -64(%rsp){%k1}   addq      $16, %rax   vmovaps   -64(%rsp), %zmm3   kmov      %k1, %r9d   popcnt    %r9d, %ecx   movl      %ecx, %ecx   shll      %cl, %r10d   xorl      $-1, %r10d   kmov      %r10d, %k2   vpackstorelps %zmm3, (%rdi,%r8,4){%k2}   vmovaps   -64(%rsp), %zmm4   nop   vpackstorehps %zmm4, 64(%rdi,%r8,4){%k2}   addq      %rcx, %r8   cmpq      $992, %rax   jb        ..B1.6