2 · Grille, blocs, warps¶

Comment découper un problème. C'est la première décision de conception d'un noyau, et celle qu'on regrette le plus tard si on la prend mal.

2.1 Les trois niveaux, et qui décide de quoi¶

Niveau	Choisi par	Contraintes
Grille (nombre de blocs)	vous	\(\le 2^{31}-1\) en `x`, \(\le 65\,535\) en `y` et `z`
Bloc (threads par bloc)	vous	\(\le 1\,024\) threads au total
Warp (32 threads)	le matériel	immuable

Le warp est invisible dans la syntaxe mais gouverne les performances. Une taille de bloc qui n'est pas un multiple de 32 gaspille des voies d'exécution : un bloc de 100 threads occupe 4 warps, dont le dernier n'a que 4 threads actifs sur 32.

Règle numéro un

La taille de bloc est toujours un multiple de 32. En pratique : 128, 256 ou 512. 256 est le défaut raisonnable.

2.2 Choisir la taille de bloc¶

Il n'y a pas de réponse universelle. Voici le raisonnement.

Trop petit (< 64) :

pas assez de warps par bloc pour cacher la latence à l'intérieur du bloc ;
on atteint vite la limite matérielle de blocs résidents par SM (32 sur Hopper), ce qui plafonne l'occupancy. 32 blocs × 32 threads = 1 024 threads, soit 50 % d'occupancy maximum.

Trop grand (1 024) :

granularité grossière : un SM ne peut héberger que 2 blocs, donc si un bloc se termine, la moitié du SM se vide ;
__syncthreads() coûte plus cher (plus de threads à attendre) ;
l'allocation de mémoire partagée par bloc devient contrainte.

La zone raisonnable : 128 à 512.

Une méthode pratique : laisser CUDA calculer l'optimum d'occupancy.

int bloc_min_grille, bloc_taille;
cudaOccupancyMaxPotentialBlockSize(&bloc_min_grille, &bloc_taille,
                                   mon_noyau, 0, 0);
// bloc_taille est la taille maximisant l'occupancy théorique

L'occupancy maximale n'est pas la performance maximale

cudaOccupancyMaxPotentialBlockSize optimise un proxy, pas votre temps d'exécution. Il donne un bon point de départ, jamais une réponse finale. Voir Performance · Occupancy.

La seule méthode fiable reste de balayer {128, 256, 512} et de mesurer.

2.3 La boucle à parcours de grille¶

Le motif le plus important de ce chapitre. Au lieu de « un thread = un élément », on écrit « chaque thread traite plusieurs éléments espacés du pas de la grille ».

__global__ void saxpy_stride(int n, float a, const float* x, float* y) {
    int pas = blockDim.x * gridDim.x;          // nombre total de threads
    for (int i = blockIdx.x * blockDim.x + threadIdx.x;
         i < n;
         i += pas) {
        y[i] = a * x[i] + y[i];
    }
}

Lancement :

int blocs = 0, threads = 256;
cudaOccupancyMaxActiveBlocksPerMultiprocessor(&blocs, saxpy_stride, threads, 0);
int nb_sm;
cudaDeviceGetAttribute(&nb_sm, cudaDevAttrMultiProcessorCount, 0);
saxpy_stride<<<blocs * nb_sm, threads>>>(n, 2.0f, d_x, d_y);

Cinq avantages, tous réels :

Taille de grille découplée de la taille du problème. On lance exactement de quoi remplir le GPU, une fois pour toutes.
Les accès restent coalescés : à chaque itération, les threads voisins touchent des adresses voisines. C'est pour cela que le pas est blockDim.x * gridDim.x et non un découpage en tranches contiguës par thread.
Réutilisation des registres entre itérations : les valeurs constantes (comme a) sont chargées une fois.
Débogage facile : on peut lancer avec <<<1, 1>>> et obtenir un programme séquentiel correct.
C'est la base des noyaux persistants, donc des megakernels (partie 8). Un noyau persistant est exactement une boucle à parcours de grille où la source de travail est une file dynamique au lieu d'un indice.

Le contre-exemple à ne pas écrire

// MAUVAIS : découpage en tranches contiguës
int debut = tid * (n / nb_threads);
for (int i = debut; i < debut + n / nb_threads; ++i) { ... }

Ici, à une itération donnée, le thread 0 lit l'élément 0 et le thread 1 lit l'élément n/nb_threads. Les accès ne sont plus coalescés du tout. C'est le découpage naturel sur CPU et le pire possible sur GPU.

2.4 Les grilles multidimensionnelles¶

dim3 permet jusqu'à 3 dimensions. C'est du sucre syntaxique — le matériel linéarise — mais qui rend certains codes bien plus lisibles.

dim3 bloc(16, 16);                                  // 256 threads
dim3 grille((largeur  + 15) / 16,
            (hauteur  + 15) / 16);
traiter_image<<<grille, bloc>>>(image, largeur, hauteur);

__global__ void traiter_image(float* img, int largeur, int hauteur) {
    int x = blockIdx.x * blockDim.x + threadIdx.x;   // colonne
    int y = blockIdx.y * blockDim.y + threadIdx.y;   // ligne
    if (x >= largeur || y >= hauteur) return;

    img[y * largeur + x] = /* ... */;                // row-major
}

Ordre de linéarisation des threads dans un bloc :

\[ \text{tid} = \text{threadIdx.z} \cdot (D_x D_y) + \text{threadIdx.y} \cdot D_x + \text{threadIdx.x} \]

Donc x varie le plus vite. Pour un bloc 16×16, le warp 0 contient les threads (x=0..15, y=0) et (x=0..15, y=1) — deux demi-lignes.

Conséquence sur la coalescence

Avec un bloc 16×16 et un tableau row-major, un warp couvre deux lignes différentes de 16 éléments. Chaque ligne fait 64 octets : deux transactions de 32 octets par ligne, donc 4 transactions au total. C'est acceptable mais pas optimal.

Un bloc 32×8 donnerait un warp couvrant exactement une ligne de 32 éléments = 128 octets contigus. Souvent 10 à 20 % plus rapide, pour un changement d'une ligne.

2.5 Le dimensionnement du travail par thread¶

Un thread doit-il traiter 1 élément ou 8 ?

Traiter plusieurs éléments par thread (thread coarsening) permet :

d'amortir le coût d'indexation et de configuration ;
de réutiliser des valeurs en registres entre éléments ;
de créer du parallélisme d'instructions (plusieurs chargements en vol simultanément), ce qui cache la latence sans augmenter l'occupancy — l'idée centrale de Volkov.

// 4 éléments par thread, avec chargement vectorisé
__global__ void saxpy4(int n4, float a, const float4* x, float4* y) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n4) {
        float4 xv = x[i];
        float4 yv = y[i];
        yv.x = a * xv.x + yv.x;
        yv.y = a * xv.y + yv.y;
        yv.z = a * xv.z + yv.z;
        yv.w = a * xv.w + yv.w;
        y[i] = yv;
    }
}

Ici on gagne deux fois : moins d'instructions et des transactions mémoire de 16 octets par thread au lieu de 4 (une instruction LDG.128 au lieu de quatre LDG.32).

Le coût : plus de registres par thread, donc potentiellement moins de warps résidents. C'est un arbitrage à mesurer. Dans une GEMM par tuiles, on va jusqu'à 8×8 = 64 éléments par thread — et c'est ce qui fait la différence entre 20 % et 90 % de cuBLAS.

2.6 Ce qui n'existe pas : la synchronisation entre blocs¶

Il faut le dire explicitement parce que tout le monde essaie.

__global__ void faux() {
    // phase 1
    calcul_partiel();
    // AUCUNE PRIMITIVE ICI ne synchronise tous les blocs
    // phase 2 : lit ce que les autres blocs ont écrit  ← indéfini
}

__syncthreads() synchronise un bloc. Point.

Les quatre contournements possibles, par ordre croissant de sophistication :

Méthode	Comment	Coût
Deux noyaux	terminer le noyau, en lancer un autre	~1,3 à 5 µs par frontière
Cooperative groups	`cooperative_groups::grid_group::sync()` avec `cudaLaunchCooperativeKernel`	exige que tous les blocs soient résidents → grille limitée
Barrière logicielle	compteur atomique + attente active	correct seulement si tous les blocs sont résidents (noyau persistant)
Megakernel	dépendances fines par compteurs, pas de barrière globale	complexité de conception

Les deux dernières lignes sont le sujet de la partie 8. Retenez pour l'instant que la frontière de noyau est la barrière globale standard, et qu'elle coûte cher.

2.7 Combien de blocs lancer, en pratique¶

Trois régimes.

Régime 1 — beaucoup plus de blocs que de SM (× 10 ou plus). C'est le cas par défaut, et il est bien : l'ordonnanceur matériel équilibre la charge tout seul, et les blocs lents sont amortis.

Régime 2 — exactement une « vague » de blocs. On lance nb_SM × blocs_par_SM blocs. C'est le régime des noyaux persistants. Bon quand chaque bloc a beaucoup de travail et qu'on veut contrôler l'ordonnancement soi-même.

Régime 3 — un peu plus d'une vague (par exemple 1,1 vague). Le pire cas. La deuxième vague ne contient que 10 % des blocs, donc 90 % du GPU est inactif pendant toute sa durée. C'est l'effet de quantification de vagues (wave quantization), et il explique des chutes de performance brutales quand la taille du problème change à peine.

148 SM, 2 blocs résidents par SM → 296 blocs par vague

296 blocs → 1,00 vague → 100 % d'efficacité
300 blocs → 1,01 vague →  51 % d'efficacité   ← catastrophe
592 blocs → 2,00 vagues → 100 % d'efficacité

Ce qu'il faut en faire

Si vous choisissez librement la taille de tuile, choisissez-la pour que le nombre de blocs soit un multiple entier du nombre de blocs par vague. cuBLAS et CUTLASS font exactement ça, et c'est une des raisons de leur avance sur les implémentations naïves.

Résumé du chapitre¶

À retenir

Taille de bloc : multiple de 32, entre 128 et 512, 256 par défaut.
La boucle à parcours de grille découple la grille du problème, préserve la coalescence, et prépare aux noyaux persistants.
En 2D, x doit indexer la dimension contiguë. Préférer 32×8 à 16×16 pour un tableau row-major.
Traiter plusieurs éléments par thread crée du parallélisme d'instructions et permet les chargements vectorisés (float4).
Il n'y a pas de synchronisation inter-blocs standard. La frontière de noyau est la barrière, et elle coûte de 1,3 à 5 µs.
Attention à la quantification de vagues : 1,01 vague coûte presque autant que 2.

Vérifiez que vous avez compris¶

Vous lancez 1 000 000 de blocs de 32 threads sur un H100. Que se passe-t-il ?

Cela fonctionne, mais c'est inefficace pour deux raisons :

32 threads par bloc = 1 warp par bloc. La limite de 32 blocs résidents par SM plafonne à 32 warps = 1 024 threads, soit 50 % d'occupancy maximum, quoi que vous fassiez.
1 million de blocs ont chacun un coût d'ordonnancement et de configuration de mémoire partagée.

Avec des blocs de 256 threads, il faudrait 125 000 blocs et l'occupancy plafonnerait à 100 %.

Pourquoi la boucle à parcours de grille utilise-t-elle blockDim.x * gridDim.x comme pas, et pas simplement blockDim.x ?

Parce que le pas doit être le nombre total de threads de la grille. Avec un pas de blockDim.x, tous les blocs traiteraient les mêmes indices : le bloc 0 et le bloc 1 feraient exactement le même travail (bug de correction, pas seulement de performance).

Le pas blockDim.x * gridDim.x garantit une partition exacte : chaque élément est traité par exactement un thread, et à chaque itération les threads consécutifs touchent des adresses consécutives.

Un noyau de convolution 2D utilise un bloc 16×16 sur une image 1920×1080. Combien de vagues sur un H100 avec 4 blocs résidents par SM ?

Blocs : \(\lceil 1920/16 \rceil \times \lceil 1080/16 \rceil = 120 \times 68 = 8\,160\).

Blocs par vague : \(132 \times 4 = 528\).

Vagues : \(8\,160 / 528 = 15{,}45\).

La fraction 0,45 de la dernière vague signifie que 55 % du GPU est inactif pendant la dernière vague, soit une perte d'efficacité globale d'environ 3,5 %. Acceptable ici parce qu'il y a beaucoup de vagues — le problème ne devient grave qu'en dessous de 3 ou 4 vagues.

Chapitre suivant : 3 · La mémoire partagée

Sources de ce chapitre¶

CUDA C++ Programming Guide — Thread Hierarchy
CUDA Pro Tip: Write Flexible Kernels with Grid-Stride Loops
CUDA C++ Best Practices Guide — Occupancy
Efficient GEMM in CUDA, CUTLASS pour la quantification de vagues et les tile schedulers.