Aller au contenu

3 · Le modèle SIMT

Trente-deux threads, une seule instruction. C'est le mécanisme d'exécution du GPU, la source de ses performances, et la source de la moitié de ses pièges.

3.1 SIMD, SIMT, et la différence

SIMD (Single Instruction, Multiple Data) est le modèle des extensions vectorielles CPU : une instruction vaddps additionne 16 float d'un coup. Le programmeur voit un vecteur.

SIMT (Single Instruction, Multiple Thread) est le modèle NVIDIA : une instruction est exécutée par 32 threads simultanément. Le programmeur voit 32 threads indépendants qui ont chacun leur threadIdx, leurs variables locales, et peuvent prendre des branches différentes.

C'est la même mécanique matérielle, exposée différemment.

SIMD (AVX) SIMT (CUDA)
Ce que le programmeur écrit _mm512_add_ps(a, b) c[i] = a[i] + b[i];
Divergence à gérer à la main par masques gérée par le matériel
Accès mémoire dispersés gather/scatter explicites naturels (mais lents)
Largeur visible explicite (16 float) invisible (mais toujours 32)

Intuition

SIMT, c'est du SIMD auquel on a ajouté un mécanisme automatique de masques d'activation. Cela rend le code beaucoup plus lisible et cache un coût que vous devrez apprendre à voir.

3.2 Le warp

Un warp est un groupe de 32 threads consécutifs, formé automatiquement à partir du threadIdx :

\[ \text{warp\_id} = \left\lfloor \frac{\text{tid}}{32} \right\rfloor, \qquad \text{lane\_id} = \text{tid} \bmod 32 \]

où tid est l'indice linéaire du thread dans son bloc.

Propriétés fondamentales :

  • le warp est l'unité d'ordonnancement. Le matériel ne planifie jamais un thread seul ;
  • tous les threads d'un warp progressent ensemble (avec les nuances de la section 3.5) ;
  • le warp est l'unité d'accès mémoire. Les 32 requêtes d'un warp sont fusionnées en un minimum de transactions ;
  • la taille de 32 est matérielle et immuable sur toutes les architectures NVIDIA depuis 2006. Chez AMD, un wavefront CDNA fait 64.

Erreur fréquente

Utiliser la variable warpSize en supposant qu'elle vaut 32 à la compilation. Sur NVIDIA c'est vrai, mais warpSize n'est pas une constante de compilation en CUDA — il faut écrire 32 explicitement si on veut qu'elle soit utilisable dans un template ou un #pragma unroll. Sur HIP, c'est warpSize = 64 sur CDNA et 32 sur RDNA, et ce piège est réel.

3.3 Le bloc et la grille

Au-dessus du warp, deux niveaux logiciels :

grille  (grid)     ← tout le lancement de noyau
 └─ bloc (block)   ← affecté à un SM, ne le quitte jamais
     └─ warp       ← 32 threads, unité d'ordonnancement
         └─ thread

Un bloc :

  • est affecté entièrement à un seul SM, pour toute sa durée de vie ;
  • peut se synchroniser en interne avec __syncthreads() ;
  • partage une allocation de mémoire partagée ;
  • a une taille choisie par le programmeur, maximum 1 024 threads.

Une grille est l'ensemble des blocs. Il n'y a aucune garantie d'ordre entre les blocs, et aucun moyen standard de les synchroniser. Deux blocs peuvent s'exécuter simultanément, ou l'un après l'autre, ou dans n'importe quel ordre.

La limite fondamentale du modèle CUDA classique

L'absence de synchronisation inter-blocs n'est pas un oubli : c'est ce qui permet à un noyau de tourner sur un GPU de 20 SM comme sur un GPU de 148 SM. Le matériel lance les blocs au fur et à mesure que des SM se libèrent.

Le prix : la seule barrière globale disponible est la fin du noyau. C'est exactement le coût que les megakernels cherchent à contourner.

Hopper a ajouté un niveau intermédiaire, le cluster : un groupe de blocs garantis co-résidents sur des SM voisins, qui peuvent se synchroniser entre eux. Voir CUDA moderne · Clusters.

3.4 La divergence

Voici le concept qui fait échouer les premiers noyaux de tout le monde.

if (threadIdx.x % 2 == 0) {
    resultat = calcul_A();   // 16 threads du warp
} else {
    resultat = calcul_B();   // les 16 autres
}

Il n'y a qu'un seul compteur ordinal par warp (nuance en 3.5). Le matériel ne peut donc pas exécuter calcul_A et calcul_B en même temps. Il fait ceci :

1. exécute calcul_A avec un masque activant les threads pairs
   → les 16 threads impairs sont présents mais leurs résultats sont jetés
2. exécute calcul_B avec le masque complémentaire
   → les 16 threads pairs sont inactifs
3. réunit les deux chemins

Coût : \(t_A + t_B\) au lieu de \(\max(t_A, t_B)\). Avec un switch à 32 branches toutes prises par un thread différent, on paie 32 fois le temps d'une branche, soit une efficacité de 1/32.

Ce qui n'est PAS de la divergence

C'est aussi important que le reste :

if (blockIdx.x % 2 == 0) { ... } else { ... }

Ici, tous les threads d'un warp prennent la même branche (ils ont le même blockIdx). Il n'y a aucune divergence, aucun coût. Le branchement est uniforme au niveau du warp.

De même :

if (i < n) { ... }   // garde de bord classique

ne diverge que dans le dernier warp du dernier bloc. Négligeable.

La règle opérationnelle

La divergence coûte quand des threads du même warp prennent des chemins différents. Rangez vos données pour que les threads voisins fassent le travail voisin. Le tri par catégorie avant traitement est un motif classique et rentable.

Un cas où la divergence est déguisée

while (donnees[i] != FIN) { ... }   // longueur variable par thread

Le warp entier tourne aussi longtemps que le thread le plus lent. Si un thread traite 1 000 éléments et les 31 autres en traitent 10, le warp coûte 1 000 itérations. C'est le problème du déséquilibre de charge intra-warp, très courant en traitement de graphes et en parcours de rayons.

3.5 Volta a changé les règles : l'exécution indépendante des threads

Jusqu'à Pascal (2016), un warp avait littéralement un compteur ordinal unique et une pile de reconvergence. Depuis Volta (2017), chaque thread a son propre compteur ordinal et sa propre pile d'appels.

Le matériel peut donc entrelacer les branches divergentes au lieu de les sérialiser strictement, et surtout : des threads divergents peuvent communiquer entre eux, ce qui était impossible avant.

La contrepartie est brutale : il n'y a plus de reconvergence implicite garantie. Ce code, valide avant Volta, est un bug depuis :

// BUG depuis Volta
if (lane < 16) {
    partage[lane] += partage[lane + 16];
}
// Rien ne garantit que l'écriture ci-dessus est visible ici
valeur = partage[0];

Il faut désormais synchroniser explicitement :

if (lane < 16) {
    partage[lane] += partage[lane + 16];
}
__syncwarp();          // ← obligatoire
valeur = partage[0];

Et toutes les primitives de warp prennent un masque explicite indiquant quels threads participent :

// Ancienne API, supprimée
int v = __shfl_down(x, 16);

// API actuelle : masque explicite
int v = __shfl_down_sync(0xffffffff, x, 16);

Limite importante

Beaucoup de code CUDA trouvé en ligne date d'avant 2017 et omet ces synchronisations. Il semble fonctionner — souvent pendant des mois — puis casse au changement de compilateur ou de carte. Si vous lisez un tutoriel de réduction sans __syncwarp() ni _sync, il est obsolète.

3.6 Les primitives de warp

Puisque les 32 threads d'un warp sont physiquement dans la même partition, ils peuvent s'échanger des registres sans passer par la mémoire. C'est extrêmement rapide (quelques cycles) et c'est la base de toutes les réductions performantes.

Primitive Ce qu'elle fait
__shfl_sync(m, v, src) chaque thread reçoit v du thread src
__shfl_up_sync(m, v, d) reçoit v du thread lane − d
__shfl_down_sync(m, v, d) reçoit v du thread lane + d
__shfl_xor_sync(m, v, k) reçoit v du thread lane XOR k
__ballot_sync(m, pred) renvoie un masque 32 bits des threads où pred est vrai
__all_sync / __any_sync réduction booléenne sur le warp
__activemask() quels threads sont actuellement actifs
__match_any_sync(m, v) quels threads ont la même valeur v
__reduce_add_sync(m, v) réduction matérielle (Ampere et plus)

Une somme sur tout un warp, en 5 instructions au lieu de 32 :

__device__ float somme_warp(float x) {
    // Réduction en arbre : 16, 8, 4, 2, 1
    for (int offset = 16; offset > 0; offset /= 2) {
        x += __shfl_down_sync(0xffffffff, x, offset);
    }
    return x;   // le résultat est dans la lane 0
}

Déroulement pour 4 threads (au lieu de 32, pour tenir sur la page) :

départ     : [a]   [b]   [c]   [d]
offset = 2 : [a+c] [b+d] [c]   [d]
offset = 1 : [a+c+b+d] …

\(\log_2(32) = 5\) étapes. Chaque étape est un échange de registre, sans aucun accès mémoire.

Pourquoi 0xffffffff ?

C'est le masque « les 32 threads participent ». Il est correct si et seulement si le warp est plein et non divergent à cet endroit. Dans un noyau avec des gardes de bord, il faut utiliser __activemask() ou, mieux, structurer le code pour que le warp soit toujours complet.

3.7 Le prédicat, ou comment le matériel évite les branchements

Pour de très courtes branches, le compilateur ne génère pas de saut du tout : il utilise la prédication. L'instruction est exécutée par tous les threads, mais son écriture est annulée pour ceux dont le prédicat est faux.

x = (cond) ? a : b;

devient en SASS quelque chose comme :

ISETP.GT P0, R2, R3      // P0 = cond
@P0  MOV R4, R5          // exécuté si P0
@!P0 MOV R4, R6          // exécuté si !P0

Aucun branchement, aucune divergence au sens strict — mais les deux chemins sont toujours payés. C'est bien pour deux instructions, catastrophique pour deux cents. Le compilateur choisit le seuil ; on peut l'influencer avec #pragma unroll ou en restructurant.

Résumé du chapitre

À retenir

  • Un warp = 32 threads (64 en wavefront AMD CDNA), unité d'ordonnancement et unité d'accès mémoire.
  • Un bloc vit sur un seul SM ; il n'existe aucune synchronisation inter-blocs standard. C'est le point de départ des megakernels.
  • La divergence coûte la somme des branches, mais seulement entre threads du même warp.
  • Depuis Volta, chaque thread a son compteur ordinal : plus de reconvergence implicite, et toutes les primitives de warp exigent un masque _sync.
  • Les échanges de registres intra-warp (__shfl_*) sont la base de toute réduction rapide.

Vérifiez que vous avez compris

Un noyau contient if (donnees[i] > 0) traiter(); où donnees est aléatoire, avec 50 % de positifs. Quelle est l'efficacité d'exécution attendue ?

Pour un warp donné, la probabilité que les 32 threads soient tous du même côté est \(2 \times 0{,}5^{32} \approx 5 \times 10^{-10}\) : essentiellement jamais. Presque tous les warps divergent, et paient les deux chemins.

Si traiter() coûte \(T\) et le chemin vide coûte 0, on paie \(T\) par warp, mais seulement la moitié des threads en profitent : efficacité ≈ 50 %.

Le remède classique : trier ou compacter les indices positifs dans un tableau séparé (stream compaction), puis lancer un noyau dense dessus. On paie un passage supplémentaire pour gagner un facteur 2 sur le traitement — rentable dès que traiter() est cher.

Pourquoi __syncthreads() dans une branche divergente est-il un comportement indéfini ?

__syncthreads() est une barrière au niveau du bloc : elle attend que tous les threads du bloc l'atteignent. Si seuls certains threads y arrivent (parce que les autres ont pris la branche else), la barrière n'est jamais satisfaite. Sur les architectures pré-Volta, cela produisait un blocage. Depuis Volta le comportement est mieux défini mais reste dangereux : la règle est toujours placer __syncthreads() dans du code atteint uniformément par tout le bloc.

Réécrivez somme_warp pour un wavefront AMD de 64 threads.

Il faut une étape de plus, et le masque n'a pas le même sens :

__device__ float somme_wave(float x) {
    for (int offset = 32; offset > 0; offset /= 2) {
        x += __shfl_down(x, offset, 64);
    }
    return x;
}

\(\log_2(64) = 6\) étapes. C'est le genre de détail qui rend un portage CUDA → HIP non trivial malgré la compatibilité syntaxique.

Chapitre suivant : 4 · La hiérarchie mémoire

Sources de ce chapitre