Aller au contenu

5 · Réductions et scans

Deux algorithmes, et presque tout le reste en dérive. La réduction est l'exercice canonique de la programmation GPU : sept versions successives, chacune plus rapide que la précédente, chacune enseignant un principe.

5.1 Ce qu'est une réduction

Combiner \(n\) éléments en un seul avec un opérateur associatif :

\[ r = x_0 \oplus x_1 \oplus \dots \oplus x_{n-1} \]

L'associativité est ce qui autorise le parallélisme : on peut regrouper comme on veut.

\[ ((x_0 \oplus x_1) \oplus (x_2 \oplus x_3)) = (x_0 \oplus (x_1 \oplus (x_2 \oplus x_3))) \]

Opérateurs usuels : somme, produit, max, min, ET/OU logique, argmax.

L'addition flottante n'est pas associative

Elle l'est mathématiquement, pas en arithmétique flottante. Une réduction parallèle ne donne donc pas le même résultat qu'une boucle séquentielle. Ce n'est en général pas un problème — l'erreur d'une réduction en arbre est même plus faible que celle d'une boucle séquentielle, en \(O(\log n)\) au lieu de \(O(n)\) — mais cela rend les résultats non reproductibles bit à bit si l'ordre varie d'une exécution à l'autre.

Voir Fondations · Nombres flottants.

Complexité. Séquentiellement : \(O(n)\) opérations, \(O(n)\) étapes. En parallèle avec \(n/2\) processeurs : \(O(n)\) opérations, \(O(\log n)\) étapes. La réduction est work-efficient : le parallélisme ne coûte aucune opération supplémentaire. C'est rare et c'est pour cela qu'elle est un si bon exemple.

5.2 Version 1 — naïve, et pourquoi elle est mauvaise

__global__ void reduce_v1(const float* entree, float* sortie) {
    extern __shared__ float sdata[];
    unsigned tid = threadIdx.x;
    unsigned i   = blockIdx.x * blockDim.x + threadIdx.x;

    sdata[tid] = entree[i];
    __syncthreads();

    for (unsigned s = 1; s < blockDim.x; s *= 2) {
        if (tid % (2 * s) == 0) {           // ← le problème
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) sortie[blockIdx.x] = sdata[0];
}

Le problème : tid % (2*s) == 0. À la première itération, seuls les threads pairs travaillent — divergence maximale dans chaque warp. Et l'opérateur modulo est lent.

À l'étape \(s\), la fraction de threads actifs est \(1/(2s)\), et les threads actifs sont dispersés. Un warp entier est mobilisé pour 1 thread utile dès \(s = 32\).

5.3 Version 2 — indexation sans divergence

for (unsigned s = 1; s < blockDim.x; s *= 2) {
    int index = 2 * s * tid;
    if (index < blockDim.x) {
        sdata[index] += sdata[index + s];
    }
    __syncthreads();
}

Maintenant les threads actifs sont les premiers, donc les warps entiers sont soit tous actifs soit tous inactifs. Plus de divergence.

Mais : sdata[index] avec index = 2*s*tid produit des conflits de banc massifs. Pour \(s = 16\), index = 32*tid → tous les threads visent le banc 0.

5.4 Version 3 — adressage séquentiel

for (unsigned s = blockDim.x / 2; s > 0; s >>= 1) {
    if (tid < s) {
        sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
}

On inverse le sens : on part de \(s = n/2\) et on descend. Les threads actifs sont 0..s-1, contigus, et ils accèdent à sdata[tid] — un banc par thread, aucun conflit.

C'est la version « correcte » de base. Elle a encore deux défauts.

5.5 Version 4 — première addition au chargement

À la première itération, la moitié des threads ne fait rien d'autre que charger. Autant leur faire additionner deux éléments tout de suite :

unsigned i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
sdata[tid] = entree[i] + entree[i + blockDim.x];
__syncthreads();

On divise par deux le nombre de blocs, donc le nombre d'étapes globales.

5.6 Version 5 — dérouler le dernier warp

Quand \(s \le 32\), il ne reste qu'un warp actif. __syncthreads() est alors inutile — mais depuis Volta, __syncwarp() reste nécessaire.

La bonne façon de l'écrire en 2026 utilise les shuffles :

__global__ void reduce_v5(const float* entree, float* sortie, int n) {
    extern __shared__ float sdata[];
    unsigned tid = threadIdx.x;
    unsigned i   = blockIdx.x * (blockDim.x * 2) + tid;

    float x = (i < n ? entree[i] : 0.0f)
            + (i + blockDim.x < n ? entree[i + blockDim.x] : 0.0f);
    sdata[tid] = x;
    __syncthreads();

    // Réduction en mémoire partagée jusqu'à 32 éléments
    for (unsigned s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) sdata[tid] += sdata[tid + s];
        __syncthreads();
    }

    // Dernier warp : shuffles, plus aucun __syncthreads()
    if (tid < 32) {
        x = sdata[tid] + sdata[tid + 32];
        #pragma unroll
        for (int off = 16; off > 0; off >>= 1)
            x += __shfl_down_sync(0xffffffff, x, off);
        if (tid == 0) sortie[blockIdx.x] = x;
    }
}

5.7 Version 6 — la version moderne, à deux niveaux

En 2026, on n'écrit plus la boucle en mémoire partagée du tout. On fait réduction de warp → partiels en mémoire partagée → réduction de warp, comme vu au chapitre précédent, le tout dans une boucle à parcours de grille :

__device__ __forceinline__ float reduce_warp(float x) {
    #pragma unroll
    for (int off = 16; off > 0; off >>= 1)
        x += __shfl_down_sync(0xffffffff, x, off);
    return x;
}

__global__ void reduce_v6(const float* entree, float* sortie, int n) {
    __shared__ float partiels[32];

    // 1. Accumulation séquentielle par thread (grid-stride)
    float somme = 0.0f;
    for (int i = blockIdx.x * blockDim.x + threadIdx.x;
         i < n;
         i += blockDim.x * gridDim.x) {
        somme += entree[i];
    }

    // 2. Réduction intra-warp
    somme = reduce_warp(somme);

    // 3. Un partiel par warp
    int lane = threadIdx.x & 31;
    int warp = threadIdx.x >> 5;
    if (lane == 0) partiels[warp] = somme;
    __syncthreads();

    // 4. Le premier warp réduit les partiels
    if (warp == 0) {
        somme = (lane < (blockDim.x >> 5)) ? partiels[lane] : 0.0f;
        somme = reduce_warp(somme);
        if (lane == 0) atomicAdd(sortie, somme);   // ou sortie[blockIdx.x]
    }
}

Ce noyau a un seul __syncthreads() et utilise atomicAdd une fois par bloc, ce qui évite un second lancement de noyau. Sur un H100 il atteint couramment 85 à 95 % de la bande passante mémoire.

L'étape 1 est la plus importante

La boucle à parcours de grille accumule séquentiellement dans un registre. Avec 132 SM × 2 048 threads = 270 336 threads et \(n = 10^9\), chaque thread accumule ~3 700 éléments avant que la moindre communication ait lieu. Le coût de la réduction en arbre devient négligeable.

C'est le principe général : maximiser le travail séquentiel par thread, minimiser les étapes de communication.

5.8 Version 7 — ne pas l'écrire

#include <cub/cub.cuh>

void reduire(const float* d_in, float* d_out, int n) {
    void*  d_temp = nullptr;
    size_t taille_temp = 0;
    cub::DeviceReduce::Sum(d_temp, taille_temp, d_in, d_out, n);
    cudaMalloc(&d_temp, taille_temp);
    cub::DeviceReduce::Sum(d_temp, taille_temp, d_in, d_out, n);
    cudaFree(d_temp);
}

CUB (partie de CCCL, livré avec le CUDA Toolkit) contient des réductions, scans, tris et sélections réglés pour chaque architecture. Il bat votre version manuelle dans la quasi-totalité des cas.

CUDA 13.1 a d'ailleurs simplifié ces APIs à une seule phase (plus besoin du double appel pour dimensionner le tampon) et ajouté des options de réduction flottante déterministe.

Pourquoi apprendre les six versions précédentes, alors ?

Parce qu'elles enseignent, dans l'ordre : la divergence, les conflits de banc, l'utilisation des threads inactifs, le déroulage, les shuffles, et la hiérarchie de communication. Ces six leçons se réappliquent à tous vos noyaux, y compris ceux pour lesquels aucune bibliothèque n'existe.

C'est aussi la structure de la lecture 9 de GPU MODE et du célèbre exposé de Mark Harris, Optimizing Parallel Reduction in CUDA, qui reste après quinze ans le meilleur document d'introduction à l'optimisation GPU.

5.9 Le scan (somme préfixe)

Le scan calcule toutes les réductions préfixes :

\[ y_i = x_0 \oplus x_1 \oplus \dots \oplus x_i \quad \text{(scan inclusif)} \]
\[ y_i = x_0 \oplus \dots \oplus x_{i-1}, \quad y_0 = \text{identité} \quad \text{(scan exclusif)} \]

C'est moins intuitif qu'une réduction et bien plus utile qu'il n'y paraît. Le scan est la primitive de :

  • l'allocation d'espace de sortie de taille variable (compaction, tri, construction de listes d'adjacence) ;
  • le tri par base (radix sort) ;
  • la construction d'histogrammes cumulés ;
  • la répartition de tokens vers les experts dans un MoE (IA · MoE) ;
  • le calcul d'offsets pour les séquences de longueur variable dans un lot (cu_seqlens de FlashAttention).

L'algorithme de Blelloch

En deux passes sur un arbre binaire, \(O(n)\) opérations, \(O(\log n)\) étapes.

Phase montante (up-sweep, identique à une réduction) :

[ 3  1  7  0  4  1  6  3 ]
[ 3  4  7  7  4  5  6  9 ]     d=0 : additionne les paires
[ 3  4  7 11  4  5  6 14 ]     d=1
[ 3  4  7 11  4  5  6 25 ]     d=2 : la racine contient la somme totale

Phase descendante (down-sweep) : on met 0 à la racine, puis à chaque niveau, l'enfant gauche reçoit la valeur du parent et l'enfant droit reçoit parent + ancienne valeur de l'enfant gauche.

[ 3  4  7 11  4  5  6  0 ]
[ 3  4  7  0  4  5  6 11 ]
[ 3  0  7  4  4 11  6 16 ]
[ 0  3  4 11 11 15 16 22 ]     ← scan exclusif

Le scan sur un warp

Avec __shfl_up_sync, en 5 étapes :

__device__ __forceinline__ float scan_inclusif_warp(float x) {
    int lane = threadIdx.x & 31;
    #pragma unroll
    for (int off = 1; off < 32; off <<= 1) {
        float v = __shfl_up_sync(0xffffffff, x, off);
        if (lane >= off) x += v;
    }
    return x;
}

Le scan sur une grille : decoupled look-back

Le scan global pose un problème que la réduction n'a pas : le bloc \(i\) a besoin de la somme de tous les blocs précédents. La solution naïve demande trois lancements de noyau.

L'algorithme decoupled look-back (Merrill & Garland, NVIDIA, 2016) le fait en un seul passage : chaque bloc publie d'abord son agrégat local, puis regarde en arrière les blocs précédents jusqu'à trouver un préfixe complet. Il atteint la bande passante mémoire — c'est-à-dire qu'un scan coûte le même prix qu'une copie.

C'est ce qu'implémente cub::DeviceScan. Ne l'écrivez pas vous-même.

cub::DeviceScan::ExclusiveSum(d_temp, taille_temp, d_in, d_out, n);

Erreur fréquente

Implémenter un scan global avec atomicAdd sur un compteur partagé, en supposant que les blocs s'exécutent dans l'ordre de blockIdx. Aucune garantie n'existe sur cet ordre. Le code fonctionne en test et casse en production. Si vous avez vraiment besoin d'un ordre, il faut le forcer avec un dynamic block index (un atomicAdd qui distribue les rangs à l'entrée du noyau), ce que fait justement le decoupled look-back.

Résumé du chapitre

À retenir

  • Une réduction parallèle : \(O(\log n)\) étapes sans opération supplémentaire.
  • Les sept versions enseignent, dans l'ordre : divergence → conflits de banc → threads inactifs → déroulage → shuffles → hiérarchie de communication → « utilisez CUB ».
  • Le motif moderne : accumulation séquentielle par thread (grid-stride), puis réduction de warp, puis partiels, puis réduction finale. Un seul __syncthreads().
  • Le scan est la primitive de l'allocation de taille variable ; il coûte le prix d'une copie mémoire avec cub::DeviceScan.
  • Ne supposez jamais un ordre d'exécution entre blocs.

Vérifiez que vous avez compris

Pourquoi la version 3 (adressage séquentiel) est-elle meilleure que la version 2, alors que les deux évitent la divergence ?

La version 2 accède à sdata[2*s*tid]. Pour \(s = 16\), l'écart entre threads consécutifs est de 32 float = 128 octets, donc tous les threads visent le même banc : conflit à 32 voies.

La version 3 accède à sdata[tid] avec tid contigu : un banc par thread, accès en un cycle. C'est un facteur 32 sur les étapes concernées.

Dans la version 6, l'atomicAdd final n'est-il pas un goulot d'étranglement ?

Non, parce qu'il n'y a qu'un atomique par bloc. Avec une grille de quelques centaines de blocs, cela fait quelques centaines d'atomiques sur la même adresse — négligeable devant les milliards d'octets lus.

Ce serait un goulot si on faisait un atomique par thread : 270 000 atomiques sérialisés sur la même ligne de cache. C'est exactement pourquoi on réduit hiérarchiquement d'abord.

Une réserve : atomicAdd flottant rend le résultat non déterministe, puisque l'ordre d'arrivée des blocs varie. Si la reproductibilité compte, écrivez sortie[blockIdx.x] et faites un second passage.

Vous devez compacter un tableau de 10⁸ éléments en gardant ceux qui vérifient un prédicat. Quelle primitive ?

Un scan exclusif sur les indicateurs (0/1) donne directement la position d'écriture de chaque élément retenu :

valeurs   : [ a  b  c  d  e  f ]
prédicat  : [ 1  0  1  1  0  1 ]
scan excl.: [ 0  1  1  2  3  3 ]   ← position de sortie
sortie    : [ a  c  d  f ]

En pratique : cub::DeviceSelect::Flagged ou cub::DeviceSelect::If, qui font exactement cela en un passage à la bande passante mémoire.

Chapitre suivant : 6 · Flux, événements et graphes

Sources de ce chapitre