cuda
Riduzione parallela (es. Come sommare una matrice)
Ricerca…
Osservazioni
L'algoritmo di riduzione parallela si riferisce in genere a un algoritmo che combina una serie di elementi, producendo un singolo risultato. I problemi tipici che rientrano in questa categoria sono:
- riassumere tutti gli elementi in un array
- trovare un massimo in un array
In generale, la riduzione parallela può essere applicata per qualsiasi operatore associativo binario, cioè (A*B)*C = A*(B*C) . Con tale operatore *, l'algoritmo di riduzione parallela raggruppa ripetutamente gli argomenti dell'array in coppie. Ogni coppia è calcolata in parallelo con altre, dimezzando la dimensione complessiva dell'array in un solo passaggio. Il processo viene ripetuto fino a quando esiste un solo elemento.
Se l'operatore è commutativo (cioè A*B = B*A ) oltre ad essere associativo, l'algoritmo può accoppiarsi in un modello diverso. Dal punto di vista teorico non fa alcuna differenza, ma in pratica fornisce un modello di accesso alla memoria migliore:
Non tutti gli operatori associativi sono commutativi, ad esempio, prendere la moltiplicazione della matrice.
Riduzione parallela a blocco singolo per operatore commutativo
L'approccio più semplice alla riduzione parallela di CUDA consiste nell'assegnare un singolo blocco per eseguire l'operazione:
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);
Questo è più fattibile quando la dimensione dei dati non è molto grande (attorno a qualche elemento di migliaia). Questo di solito accade quando la riduzione fa parte di un programma CUDA più grande. Se l'input corrisponde a blockSize sin dall'inizio, il primo ciclo for può essere completamente rimosso.
Nota che, in primo luogo, quando ci sono più elementi dei thread, aggiungiamo le cose in modo completamente indipendente. Solo quando il problema viene ridotto a blockSize , i trigger di riduzione paralleli effettivi. Lo stesso codice può essere applicato a qualsiasi altro operatore commutativo, associativo, come la moltiplicazione, il minimo, il massimo, ecc.
Si noti che l'algoritmo può essere reso più veloce, ad esempio utilizzando una riduzione parallela a livello di curvatura.
Riduzione parallela a blocco singolo per operatore non commutativo
La riduzione parallela per un operatore non commutativo è un po 'più complicata, rispetto alla versione commutativa. Nell'esempio usiamo ancora un'aggiunta su interi per la semplicità. Potrebbe essere sostituito, ad esempio, con la moltiplicazione della matrice che è realmente non commutativa. Nota, nel fare ciò, 0 dovrebbe essere sostituito da un elemento neutro della moltiplicazione, cioè una matrice di 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);
Nel primo ciclo while viene eseguito finché ci sono più elementi di input dei thread. In ogni iterazione viene eseguita una singola riduzione e il risultato viene compresso nella prima metà dell'array shArr . La seconda metà viene quindi riempita con nuovi dati.
Una volta caricati tutti i dati da gArr , viene eseguito il secondo ciclo. Ora non comprimiamo più il risultato (che costa un __syncthreads() ) extra. In ogni passo il thread n accede al 2*n -esimo elemento attivo e lo aggiunge con 2*n+1 -esimo elemento:
Esistono molti modi per ottimizzare ulteriormente questo semplice esempio, ad esempio attraverso la riduzione del livello di curvatura e rimuovendo i conflitti di memoria del banco condiviso.
Riduzione parallela multiblocco per operatore commutativo
L'approccio a blocchi multipli per la riduzione parallela di CUDA rappresenta una sfida aggiuntiva rispetto all'approccio a blocco singolo, poiché i blocchi sono limitati nella comunicazione. L'idea è di lasciare che ogni blocco calcoli una parte dell'array di input e quindi un blocco finale per unire tutti i risultati parziali. Per fare ciò è possibile avviare due kernel, creando implicitamente un punto di sincronizzazione a livello di griglia.
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;
}
Uno idealmente vuole lanciare abbastanza blocchi per saturare tutti i multiprocessori sulla GPU in piena occupazione. Il superamento di questo numero, in particolare il lancio di tutti i thread quanti sono gli elementi nell'array, è controproducente. In questo modo non aumenta più la potenza di calcolo raw, ma impedisce l'utilizzo del primo loop molto efficiente.
È anche possibile ottenere lo stesso risultato usando un singolo kernel, con l'aiuto di last-block guard :
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;
}
Si noti che il kernel può essere reso più veloce, ad esempio usando una riduzione parallela a livello di curvatura.
Riduzione parallela multiblocco per operatore non commutabile
L'approccio a blocchi multipli per la riduzione parallela è molto simile all'approccio a blocco singolo. L'array di input globale deve essere suddiviso in sezioni, ciascuna ridotta da un singolo blocco. Quando si ottiene un risultato parziale da ciascun blocco, un blocco finale li riduce per ottenere il risultato finale.
-
sumNoncommSingleBlockè spiegato più dettagliatamente nell'esempio di riduzione a blocco singolo. -
lastBlockaccetta solo l'ultimo blocco che lo raggiunge. Se vuoi evitare questo, puoi dividere il kernel in due chiamate separate.
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);
}
Uno idealmente vuole lanciare abbastanza blocchi per saturare tutti i multiprocessori sulla GPU in piena occupazione. Il superamento di questo numero, in particolare il lancio di tutti i thread quanti sono gli elementi nell'array, è controproducente. In questo modo non aumenta più la potenza di calcolo raw, ma impedisce l'utilizzo del primo loop molto efficiente.
Riduzione parallela a singolo ordito per operatore commutativo
A volte la riduzione deve essere eseguita su una scala molto piccola, come parte di un kernel CUDA più grande. Supponiamo ad esempio che i dati di input contengano esattamente 32 elementi: il numero di fili in un ordito. In tale scenario può essere assegnato un singolo ordito per eseguire la riduzione. Dato che il warp viene eseguito in una perfetta sincronizzazione, molte istruzioni __syncthreads() possono essere rimosse - se confrontate con una riduzione a livello di blocco.
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 è preferibilmente un array nella memoria condivisa. Il valore dovrebbe essere lo stesso per tutti i thread nel warp. Se sumCommSingleWarp è chiamato da più shArr , shArr dovrebbe essere diverso tra i warp (lo stesso all'interno di ogni warp).
L'argomento shArr è contrassegnato come volatile per garantire che le operazioni sull'array vengano effettivamente eseguite dove indicato. Altrimenti, l'assegnazione ripetitiva a shArr[idx] può essere ottimizzata come assegnazione a un registro, poiché solo l'asserzione finale è un vero e proprio negozio di shArr . Quando ciò accade, i compiti immediati non sono visibili ad altri thread, producendo risultati errati. Nota che puoi passare un normale array non volatile come argomento di uno volatile, come quando passi non-const come parametro const.
Se non ci si preoccupa del contenuto di shArr[1..31] dopo la riduzione, si può semplificare ulteriormente il codice:
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];
}
In questa configurazione abbiamo rimosso molti if le condizioni. I thread aggiuntivi eseguono alcune aggiunte non necessarie, ma non ci interessa più i contenuti che producono. Poiché gli orditi vengono eseguiti in modalità SIMD, in realtà non risparmiamo in tempo facendo in modo che quei thread non facciano nulla. D'altra parte, la valutazione delle condizioni richiede tempi relativamente lunghi, poiché il corpo di queste affermazioni if è così piccolo. L'istruzione if iniziale può essere rimossa anche se shArr[32..47] è riempito con 0.
La riduzione del livello di curvatura può essere utilizzata anche per aumentare la riduzione a livello di blocco:
__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'argomento &r[idx & ~(warpSize-1)] è fondamentalmente r + warpIdx*32 . Questo suddivide efficacemente l'array r in blocchi di 32 elementi, e ciascun blocco viene assegnato a distorsione separata.
Riduzione parallela a singolo ordito per operatore non commutabile
A volte la riduzione deve essere eseguita su una scala molto piccola, come parte di un kernel CUDA più grande. Supponiamo ad esempio che i dati di input contengano esattamente 32 elementi: il numero di fili in un ordito. In tale scenario può essere assegnato un singolo ordito per eseguire la riduzione. Dato che il warp viene eseguito in una perfetta sincronizzazione, molte istruzioni __syncthreads() possono essere rimosse - se confrontate con una riduzione a livello di blocco.
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];
}
shArr è preferibilmente un array nella memoria condivisa. Il valore dovrebbe essere lo stesso per tutti i thread nel warp. Se sumCommSingleWarp è chiamato da più shArr , shArr dovrebbe essere diverso tra i warp (lo stesso all'interno di ogni warp).
L'argomento shArr è contrassegnato come volatile per garantire che le operazioni sull'array vengano effettivamente eseguite dove indicato. Altrimenti, l'assegnazione ripetitiva a shArr[idx] può essere ottimizzata come assegnazione a un registro, poiché solo l'asserzione finale è un vero e proprio negozio di shArr . Quando ciò accade, i compiti immediati non sono visibili ad altri thread, producendo risultati errati. Nota che puoi passare un normale array non volatile come argomento di uno volatile, come quando passi non-const come parametro const.
Se uno non si preoccupa del contenuto finale di shArr[1..31] e può pad shArr[32..47] con zeri, si può semplificare il codice di cui sopra:
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];
}
In questa configurazione abbiamo rimosso tutto if le condizioni, che costituiscono circa la metà delle istruzioni. I thread aggiuntivi eseguono alcune aggiunte non necessarie, memorizzando il risultato in celle di shArr che alla fine non hanno alcun impatto sul risultato finale. Poiché gli orditi vengono eseguiti in modalità SIMD, in realtà non risparmiamo in tempo facendo in modo che quei thread non facciano nulla.
Riduzione parallela a singolo ordito utilizzando solo registri
In genere, la riduzione viene eseguita su un array globale o condiviso. Tuttavia, quando la riduzione viene eseguita su una scala molto piccola, come parte di un kernel CUDA più grande, può essere eseguita con un singolo warp. Quando ciò accade, su Keppler o architetture superiori (CC> = 3.0), è possibile utilizzare le funzioni warp-shuffle per evitare di utilizzare la memoria condivisa.
Supponiamo ad esempio che ogni filo in un curvino abbia un unico valore di dati di input. Tutti i thread hanno 32 elementi, che dobbiamo riassumere (o eseguire altre operazioni associative)
__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);
}
Questa versione funziona sia per operatori commutativi che non commutativi.