cuda => Réduction parallèle (par exemple comment additionner un tableau)

Remarques

L'algorithme de réduction parallèle fait généralement référence à un algorithme qui combine un ensemble d'éléments, produisant un seul résultat. Les problèmes typiques qui entrent dans cette catégorie sont les suivants:

résumer tous les éléments d'un tableau
trouver un maximum dans un tableau

En général, la réduction parallèle peut être appliquée à tout opérateur associatif binaire, à savoir (A*B)*C = A*(B*C) . Avec un tel opérateur *, l'algorithme de réduction parallèle regroupe à répétition les arguments du tableau par paires. Chaque paire est calculée en parallèle avec d'autres, réduisant de moitié la taille globale de la matrice en une seule étape. Le processus est répété jusqu'à ce qu'un seul élément existe.

Si l'opérateur est commutatif (c.-à-d. A*B = B*A ) en plus d'être associatif, l'algorithme peut se coupler selon un modèle différent. Du point de vue théorique, cela ne fait aucune différence, mais en pratique, cela donne un meilleur modèle d'accès à la mémoire:

Tous les opérateurs associatifs ne sont pas commutatifs - prenons par exemple la multiplication matricielle.

Réduction parallèle monobloc pour opérateur commutatif

L'approche la plus simple pour la réduction parallèle dans CUDA consiste à affecter un seul bloc pour effectuer la tâche:

static const int arraySize = 10000;
static const int blockSize = 1024;

__global__ void sumCommSingleBlock(const int *a, int *out) {
    int idx = threadIdx.x;
    int sum = 0;
    for (int i = idx; i < arraySize; i += blockSize)
        sum += a[i];
    __shared__ int r[blockSize];
    r[idx] = sum;
    __syncthreads();
    for (int size = blockSize/2; size>0; size/=2) { //uniform
        if (idx<size)
            r[idx] += r[idx+size];
        __syncthreads();
    }
    if (idx == 0)
        *out = r[0];
}

...

sumCommSingleBlock<<<1, blockSize>>>(dev_a, dev_out);

Cela est plus facile lorsque la taille des données n'est pas très grande (autour de quelques milliers d'éléments). Cela se produit généralement lorsque la réduction fait partie d'un programme CUDA plus important. Si l'entrée correspond à blockSize dès le début, la première for la boucle peut être complètement enlevée.

Notez que dans la première étape, quand il y a plus d'éléments que de threads, nous ajoutons les choses complètement indépendamment. Ce n'est que lorsque le problème est réduit à blockSize la réduction parallèle réelle se déclenche. Le même code peut être appliqué à tout autre opérateur associatif commutatif, tel que multiplication, minimum, maximum, etc.

Notez que l'algorithme peut être réalisé plus rapidement, par exemple en utilisant une réduction parallèle au niveau de la chaîne.

Réduction parallèle monobloc pour opérateur non commutatif

Faire de la réduction parallèle pour un opérateur non-commutatif est un peu plus complexe que la version commutative. Dans l'exemple, nous utilisons encore un ajout sur les entiers pour la simplicité. Il pourrait être remplacé, par exemple, par la multiplication matricielle qui est vraiment non-commutative. Notez que ce faisant, 0 doit être remplacé par un élément neutre de la multiplication, c’est-à-dire une matrice d’identité.

static const int arraySize = 1000000;
static const int blockSize = 1024;

__global__ void sumNoncommSingleBlock(const int *gArr, int *out) {
    int thIdx = threadIdx.x;
    __shared__ int shArr[blockSize*2];
    __shared__ int offset;
    shArr[thIdx] = thIdx<arraySize ? gArr[thIdx] : 0;
    if (thIdx == 0)
        offset = blockSize;
    __syncthreads();
    while (offset < arraySize) { //uniform
        shArr[thIdx + blockSize] = thIdx+offset<arraySize ? gArr[thIdx+offset] : 0;
        __syncthreads();
        if (thIdx == 0)
            offset += blockSize;
        int sum = shArr[2*thIdx] + shArr[2*thIdx+1];
        __syncthreads();
        shArr[thIdx] = sum;
    }
    __syncthreads();
    for (int stride = 1; stride<blockSize; stride*=2) { //uniform
        int arrIdx = thIdx*stride*2;
        if (arrIdx+stride<blockSize)
            shArr[arrIdx] += shArr[arrIdx+stride];
        __syncthreads();
    }   
    if (thIdx == 0)
        *out = shArr[0];
}

...

sumNoncommSingleBlock<<<1, blockSize>>>(dev_a, dev_out);

Dans la première, la boucle while s'exécute tant qu'il y a plus d'éléments d'entrée que de threads. Dans chaque itération, une seule réduction est effectuée et le résultat est compressé dans la première moitié du tableau shArr . La seconde moitié est ensuite remplie de nouvelles données.

Une fois toutes les données chargées depuis gArr , la deuxième boucle est gArr . Maintenant, nous ne compressons plus le résultat (ce qui coûte un __syncthreads() supplémentaire). A chaque étape, le thread n accède au 2*n ème élément actif et l'ajoute avec 2*n+1 -ème élément:

Il existe de nombreuses manières d’optimiser cet exemple simple, par exemple en réduisant le niveau de distorsion et en supprimant les conflits de banque de mémoire partagée.

Réduction parallèle multiblocs pour opérateur commutatif

Une approche multi-blocs pour la réduction parallèle de CUDA pose un défi supplémentaire, comparé à l'approche à un seul bloc, car les blocs sont limités dans la communication. L'idée est de laisser chaque bloc calculer une partie du tableau d'entrée, puis de disposer d'un dernier bloc pour fusionner tous les résultats partiels. Pour ce faire, il est possible de lancer deux noyaux, créant implicitement un point de synchronisation à l’échelle de la grille.

static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24; //this number is hardware-dependent; usually #SM*2 is a good number.

__global__ void sumCommMultiBlock(const int *gArr, int arraySize, int *gOut) {
    int thIdx = threadIdx.x;
    int gthIdx = thIdx + blockIdx.x*blockSize;
    const int gridSize = blockSize*gridDim.x;
    int sum = 0;
    for (int i = gthIdx; i < arraySize; i += gridSize)
        sum += gArr[i];
    __shared__ int shArr[blockSize];
    shArr[thIdx] = sum;
    __syncthreads();
    for (int size = blockSize/2; size>0; size/=2) { //uniform
        if (thIdx<size)
            shArr[thIdx] += shArr[thIdx+size];
        __syncthreads();
    }
    if (thIdx == 0)
        gOut[blockIdx.x] = shArr[0];
}

__host__ int sumArray(int* arr) {
    int* dev_arr;
    cudaMalloc((void**)&dev_arr, wholeArraySize * sizeof(int));
    cudaMemcpy(dev_arr, arr, wholeArraySize * sizeof(int), cudaMemcpyHostToDevice);

    int out;
    int* dev_out;
    cudaMalloc((void**)&dev_out, sizeof(int)*gridSize);
    
    sumCommMultiBlock<<<gridSize, blockSize>>>(dev_arr, wholeArraySize, dev_out);
    //dev_out now holds the partial result
    sumCommMultiBlock<<<1, blockSize>>>(dev_out, gridSize, dev_out);
    //dev_out[0] now holds the final result
    cudaDeviceSynchronize();
    
    cudaMemcpy(&out, dev_out, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(dev_arr);
    cudaFree(dev_out);
    return out;
}

On souhaite idéalement lancer suffisamment de blocs pour saturer tous les multiprocesseurs du processeur graphique à pleine capacité. Dépasser ce nombre - en particulier, lancer autant de threads qu'il y a d'éléments dans le tableau - est contre-productif. Cela n'augmente plus la puissance de calcul brute, mais empêche d'utiliser la première boucle très efficace.

Il est également possible d’obtenir le même résultat en utilisant un seul noyau, à l’aide de la protection du dernier bloc :

static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24;

__device__ bool lastBlock(int* counter) {
    __threadfence(); //ensure that partial result is visible by all blocks
    int last = 0;
    if (threadIdx.x == 0)
        last = atomicAdd(counter, 1);
    return __syncthreads_or(last == gridDim.x-1);
}    

__global__ void sumCommMultiBlock(const int *gArr, int arraySize, int *gOut, int* lastBlockCounter) {
    int thIdx = threadIdx.x;
    int gthIdx = thIdx + blockIdx.x*blockSize;
    const int gridSize = blockSize*gridDim.x;
    int sum = 0;
    for (int i = gthIdx; i < arraySize; i += gridSize)
        sum += gArr[i];
    __shared__ int shArr[blockSize];
    shArr[thIdx] = sum;
    __syncthreads();
    for (int size = blockSize/2; size>0; size/=2) { //uniform
        if (thIdx<size)
            shArr[thIdx] += shArr[thIdx+size];
        __syncthreads();
    }
    if (thIdx == 0)
        gOut[blockIdx.x] = shArr[0];
    if (lastBlock(lastBlockCounter)) {
        shArr[thIdx] = thIdx<gridSize ? gOut[thIdx] : 0;
        __syncthreads();
        for (int size = blockSize/2; size>0; size/=2) { //uniform
            if (thIdx<size)
                shArr[thIdx] += shArr[thIdx+size];
            __syncthreads();
        }
        if (thIdx == 0)
            gOut[0] = shArr[0];            
    }
}

__host__ int sumArray(int* arr) {
    int* dev_arr;
    cudaMalloc((void**)&dev_arr, wholeArraySize * sizeof(int));
    cudaMemcpy(dev_arr, arr, wholeArraySize * sizeof(int), cudaMemcpyHostToDevice);

    int out;
    int* dev_out;
    cudaMalloc((void**)&dev_out, sizeof(int)*gridSize);
    
    int* dev_lastBlockCounter;
    cudaMalloc((void**)&dev_lastBlockCounter, sizeof(int));
    cudaMemset(dev_lastBlockCounter, 0, sizeof(int));
    
    sumCommMultiBlock<<<gridSize, blockSize>>>(dev_arr, wholeArraySize, dev_out, dev_lastBlockCounter);
    cudaDeviceSynchronize();
    
    cudaMemcpy(&out, dev_out, sizeof(int), cudaMemcpyDeviceToHost);
    cudaFree(dev_arr);
    cudaFree(dev_out);
    return out;
}

Notez que le noyau peut être rendu plus rapide, par exemple en utilisant une réduction parallèle au niveau de la chaîne.

Réduction parallèle multibloc pour un opérateur non commutatif

L'approche à plusieurs blocs pour la réduction parallèle est très similaire à l'approche à un seul bloc. Le tableau d'entrée global doit être divisé en sections, chacune étant réduite d'un seul bloc. Lorsqu'un résultat partiel de chaque bloc est obtenu, un bloc final les réduit pour obtenir le résultat final.

sumNoncommSingleBlock est expliqué plus en détail dans l'exemple de réduction de bloc unique.
lastBlock n'accepte que le dernier bloc atteint. Si vous voulez éviter cela, vous pouvez diviser le noyau en deux appels distincts.

static const int wholeArraySize = 100000000;
static const int blockSize = 1024;
static const int gridSize = 24; //this number is hardware-dependent; usually #SM*2 is a good number.

__device__ bool lastBlock(int* counter) {
    __threadfence(); //ensure that partial result is visible by all blocks
    int last = 0;
    if (threadIdx.x == 0)
        last = atomicAdd(counter, 1);
    return __syncthreads_or(last == gridDim.x-1);
}

__device__ void sumNoncommSingleBlock(const int* gArr, int arraySize, int* out) {
    int thIdx = threadIdx.x;
    __shared__ int shArr[blockSize*2];
    __shared__ int offset;
    shArr[thIdx] = thIdx<arraySize ? gArr[thIdx] : 0;
    if (thIdx == 0)
        offset = blockSize;
    __syncthreads();
    while (offset < arraySize) { //uniform
        shArr[thIdx + blockSize] = thIdx+offset<arraySize ? gArr[thIdx+offset] : 0;
        __syncthreads();
        if (thIdx == 0)
            offset += blockSize;
        int sum = shArr[2*thIdx] + shArr[2*thIdx+1];
        __syncthreads();
        shArr[thIdx] = sum;
    }
    __syncthreads();
    for (int stride = 1; stride<blockSize; stride*=2) { //uniform
        int arrIdx = thIdx*stride*2;
        if (arrIdx+stride<blockSize)
            shArr[arrIdx] += shArr[arrIdx+stride];
        __syncthreads();
    }   
    if (thIdx == 0)
        *out = shArr[0];
}

__global__ void sumNoncommMultiBlock(const int* gArr, int* out, int* lastBlockCounter) {
    int arraySizePerBlock = wholeArraySize/gridSize;
    const int* gArrForBlock = gArr+blockIdx.x*arraySizePerBlock;
    int arraySize = arraySizePerBlock;
    if (blockIdx.x == gridSize-1)
        arraySize = wholeArraySize - blockIdx.x*arraySizePerBlock;
    sumNoncommSingleBlock(gArrForBlock, arraySize, &out[blockIdx.x]);
    if (lastBlock(lastBlockCounter))
        sumNoncommSingleBlock(out, gridSize, out);    
}

Réduction parallèle à une seule chaîne pour un opérateur commutatif

Parfois, la réduction doit être effectuée à très petite échelle, dans le cadre d'un noyau CUDA plus important. Supposons, par exemple, que les données d’entrée contiennent exactement 32 éléments - le nombre de threads dans une chaîne. Dans un tel scénario, une seule chaîne peut être affectée à la réduction. Étant donné que warp s'exécute dans une synchronisation parfaite, de nombreuses instructions __syncthreads() peuvent être supprimées - par rapport à une réduction au niveau du bloc.

static const int warpSize = 32;

__device__ int sumCommSingleWarp(volatile int* shArr) {
    int idx = threadIdx.x % warpSize; //the lane index in the warp
    if (idx<16) shArr[idx] += shArr[idx+16];
    if (idx<8) shArr[idx] += shArr[idx+8];
    if (idx<4) shArr[idx] += shArr[idx+4];
    if (idx<2) shArr[idx] += shArr[idx+2];
    if (idx==0) shArr[idx] += shArr[idx+1];
    return shArr[0];
}

shArr est de préférence un tableau en mémoire partagée. La valeur doit être la même pour tous les threads du Warp. Si sumCommSingleWarp est appelé par plusieurs shArr , shArr devrait être différent entre les warps (même dans chaque warp).

L'argument shArr est marqué comme volatile pour garantir que les opérations sur le tableau sont réellement effectuées là où elles sont indiquées. Sinon, l'affectation répétitive à shArr[idx] peut être optimisée en tant qu'affectation à un registre, seule l'assignation finale étant un magasin réel pour shArr . Lorsque cela se produit, les affectations immédiates ne sont pas visibles pour les autres threads, ce qui produit des résultats incorrects. Notez que vous pouvez passer un tableau non-volatile normal en tant qu'argument de volatile, comme lorsque vous transmettez un paramètre non-const comme paramètre const.

Si l'on ne se soucie pas du contenu de shArr[1..31] après la réduction, on peut simplifier davantage le code:

static const int warpSize = 32;

__device__ int sumCommSingleWarp(volatile int* shArr) {
    int idx = threadIdx.x % warpSize; //the lane index in the warp
    if (idx<16) {
      shArr[idx] += shArr[idx+16];
      shArr[idx] += shArr[idx+8];
      shArr[idx] += shArr[idx+4];
      shArr[idx] += shArr[idx+2];
      shArr[idx] += shArr[idx+1];
    }
    return shArr[0];
}

Dans cette configuration , nous avons supprimé beaucoup if les conditions. Les threads supplémentaires effectuent des ajouts inutiles, mais nous ne nous soucions plus du contenu qu'ils produisent. Comme les warps s'exécutent en mode SIMD, nous ne réalisons pas vraiment de gain de temps en ne faisant rien. D'un autre côté, l'évaluation des conditions prend beaucoup de temps, puisque le corps de celles if ci est si petit. L'instruction if initiale peut également être supprimée si shArr[32..47] est shArr[32..47] par 0.

La réduction au niveau de la chaîne peut également être utilisée pour augmenter la réduction au niveau du bloc:

__global__ void sumCommSingleBlockWithWarps(const int *a, int *out) {
    int idx = threadIdx.x;
    int sum = 0;
    for (int i = idx; i < arraySize; i += blockSize)
        sum += a[i];
    __shared__ int r[blockSize];
    r[idx] = sum;
    sumCommSingleWarp(&r[idx & ~(warpSize-1)]);
    __syncthreads();
    if (idx<warpSize) { //first warp only
        r[idx] = idx*warpSize<blockSize ? r[idx*warpSize] : 0;
        sumCommSingleWarp(r);
        if (idx == 0)
            *out = r[0];
    }
}

L'argument &r[idx & ~(warpSize-1)] est essentiellement r + warpIdx*32 . Cela divise effectivement le tableau r en morceaux de 32 éléments, et chaque segment est assigné à séparer la chaîne.

Réduction parallèle à une seule chaîne pour opérateur non commutatif

static const int warpSize = 32;

__device__ int sumNoncommSingleWarp(volatile int* shArr) {
    int idx = threadIdx.x % warpSize; //the lane index in the warp
    if (idx%2 == 0) shArr[idx] += shArr[idx+1];
    if (idx%4 == 0) shArr[idx] += shArr[idx+2];
    if (idx%8 == 0) shArr[idx] += shArr[idx+4];
    if (idx%16 == 0) shArr[idx] += shArr[idx+8];
    if (idx == 0) shArr[idx] += shArr[idx+16];
    return shArr[0];
}

Si l'on ne se soucie pas du contenu final de shArr[1..31] et que l'on peut shArr[32..47] avec des zéros, on peut simplifier le code ci-dessus:

static const int warpSize = 32;

__device__ int sumNoncommSingleWarpPadded(volatile int* shArr) {
    //shArr[32..47] == 0
    int idx = threadIdx.x % warpSize; //the lane index in the warp
    shArr[idx] += shArr[idx+1];
    shArr[idx] += shArr[idx+2];
    shArr[idx] += shArr[idx+4];
    shArr[idx] += shArr[idx+8];
    shArr[idx] += shArr[idx+16];
    return shArr[0];
}

Dans cette configuration , nous avons supprimé toutes if les conditions, qui constituent environ la moitié des instructions. Les threads supplémentaires effectuent des ajouts inutiles, stockant le résultat dans des cellules de shArr qui n'ont finalement aucun impact sur le résultat final. Comme les warps s'exécutent en mode SIMD, nous ne réalisons pas vraiment de gain de temps en ne faisant rien.

Réduction parallèle à une seule chaîne à l'aide de registres uniquement

En règle générale, la réduction est effectuée sur un tableau global ou partagé. Cependant, lorsque la réduction est effectuée à très petite échelle, dans le cadre d’un noyau CUDA plus important, elle peut être effectuée avec une seule chaîne. Lorsque cela se produit, sur les architectures Keppler ou supérieures (CC> = 3.0), il est possible d’utiliser des fonctions Warp-Shuffle pour éviter d’utiliser la mémoire partagée.

Supposons, par exemple, que chaque thread dans une chaîne contient une seule valeur de données d'entrée. Tous les threads ensemble ont 32 éléments, que nous devons résumer (ou effectuer d'autres opérations associatives)

__device__ int sumSingleWarpReg(int value) {
    value += __shfl_down(value, 1);
    value += __shfl_down(value, 2);
    value += __shfl_down(value, 4);
    value += __shfl_down(value, 8);
    value += __shfl_down(value, 16);
    return __shfl(value,0);
}

Cette version fonctionne à la fois pour les opérateurs commutatifs et non commutatifs.

Modified text is an extract of the original Stack Overflow Documentation

Sous licence CC BY-SA 3.0

Non affilié à Stack Overflow

cuda
Réduction parallèle (par exemple comment additionner un tableau)

Recherche…