cuda
Parallell reduktion (t.ex. hur man summerar en matris)
Sök…
Anmärkningar
Parallell reduktionsalgoritm avser vanligtvis en algoritm som kombinerar en rad element, vilket ger ett enda resultat. Typiska problem som faller inom denna kategori är:
- sammanfattar alla element i en matris
- hitta ett maximum i en matris
I allmänhet kan den parallella reduktionen tillämpas för alla binära associerande operatörer , dvs (A*B)*C = A*(B*C) . Med en sådan operatör * grupperar den parallella reduktionsalgoritmen repeterat matrisargument i par. Varje par beräknas parallellt med andra, och halverar den totala matrisstorleken i ett steg. Processen upprepas tills endast ett enda element finns.
Om operatören är kommutativ (dvs. A*B = B*A ) förutom att den är associerande kan algoritmen kopplas i ett annat mönster. Ur teoretisk synvinkel gör det ingen skillnad, men i praktiken ger det ett bättre minneåtkomstmönster:
Inte alla associativa operatörer är kommutativa - ta matrismultiplikation till exempel.
Parallell reduktion med ett block för kommutativ operatör
Den enklaste metoden för parallell reduktion i CUDA är att tilldela ett enda block för att utföra uppgiften:
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);
Detta är mest genomförbart när datastorleken inte är särskilt stor (cirka några tusen element). Detta händer vanligtvis när reduktionen är en del av något större CUDA-program. Om ingången matchar blockSize från början kan den första for slingan tas bort helt.
Observera att i det första steget, när det finns fler element än trådar, lägger vi till saker helt oberoende. Först när problemet reduceras till blockSize den faktiska parallella reduktionen. Samma kod kan tillämpas på alla andra kommutativa, associerande operatörer, såsom multiplikation, minimum, max osv.
Observera att algoritmen kan göras snabbare, till exempel genom att använda en parallell reduktion av varpnivån.
Parallell reduktion med ett block för icke-kommutativ operatör
Att göra parallellminskning för en icke-kommutativ operatör är lite mer involverad jämfört med den kommutativa versionen. I exemplet använder vi fortfarande ett tillägg över heltal för enkelhets skull. Det kan till exempel ersättas med matrismultiplikation som verkligen är icke-kommutativ. Observera att när du gör det bör 0 ersättas med ett neutralt element i multiplikationen - dvs. en identitetsmatris.
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);
I det första medan loopen körs så länge det finns fler ingångselement än trådar. I varje iteration utförs en enda reduktion och resultatet komprimeras till den första halvan av shArr matrisen. Den andra halvan fylls sedan med ny data.
När alla data laddas från gArr den andra slingan. Nu komprimerar vi inte längre resultatet (vilket kostar en extra __syncthreads() ). I varje steg får tråden n åtkomst till det 2*n aktiva elementet och lägger till det med 2*n+1 :
Det finns många sätt att ytterligare optimera detta enkla exempel, t.ex. genom att minska varpnivån och genom att ta bort bankkonflikter med delat minne.
Parallell reduktion med flera block för kommutativ operatör
Multi-block-strategi för parallell reduktion av CUDA utgör en extra utmaning jämfört med en-block-strategi, eftersom block är begränsade i kommunikation. Tanken är att låta varje block beräkna en del av inmatningsfältet och sedan ha ett sista block för att slå samman alla delresultat. För att göra det kan man starta två kärnor genom att implicit skapa en rutnätbredd synkroniseringspunkt.
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;
}
Man vill helst starta tillräckligt med block för att mätta alla multiprocessorer på GPU vid full beläggning. Överskridande av detta antal - särskilt att starta så många trådar som det finns element i matrisen - är motproduktivt. Om du gör det ökar den råa datorkraften inte längre utan förhindrar att du använder den mycket effektiva första slingan.
Det är också möjligt att uppnå samma resultat med en enda kärna med hjälp av sista blockblock :
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;
}
Observera att kärnan kan göras snabbare, till exempel genom att använda en parallellreduktion av varpnivån.
Parallell reduktion med flera block för icke-kommutativ operatör
Multi-block-strategi för parallell reduktion liknar mycket ett block-tillvägagångssätt. Den globala inmatningsgruppen måste delas upp i sektioner, var och en minskas med ett enda block. När ett partiellt resultat från varje block erhålls, reducerar ett slutligt block dessa för att erhålla det slutliga resultatet.
-
sumNoncommSingleBlockförklaras mer i detalj i exemplet med enstaka blockering. -
lastBlockaccepterar bara det sista blocket som når det. Om du vill undvika detta kan du dela kärnan i två separata samtal.
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);
}
Man vill helst starta tillräckligt med block för att mätta alla multiprocessorer på GPU vid full beläggning. Överskridande av detta antal - särskilt att starta så många trådar som det finns element i matrisen - är motproduktivt. Om du gör det ökar den råa datorkraften inte längre utan förhindrar att du använder den mycket effektiva första slingan.
Parallell reduktion med en varp för kommutativ operatör
Ibland måste reduktionen utföras i mycket liten skala, som en del av en större CUDA-kärna. Anta till exempel att inmatningsdata har exakt 32 element - antalet trådar i en varp. I ett sådant scenario kan ett enda varp tilldelas för att utföra reduktionen. Med tanke på att varp körs i en perfekt synk kan många __syncthreads() instruktioner tas bort - jämfört med en blocknivåreduktion.
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 är företrädesvis en matris i delat minne. Värdet ska vara detsamma för alla trådar i varpen. Om sumCommSingleWarp kallas av flera varp ska shArr vara olika mellan varp (samma inom varje varp).
Argumentet shArr markeras som volatile att säkerställa att operationer på arrayen faktiskt utförs där det anges. Annars kan den repetitiva tilldelningen till shArr[idx] optimeras som en tilldelning till ett register, där endast slutlig tilldelning är en verklig butik för shArr . När det händer är de omedelbara tilldelningarna inte synliga för andra trådar, vilket ger felaktiga resultat. Observera att du kan passera en normal icke-flyktig matris som ett argument för flyktig, samma som när du passerar non-const som en const-parameter.
Om man inte bryr sig om innehållet i shArr[1..31] efter minskningen kan man förenkla koden ytterligare:
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];
}
I den här installationen tog vi bort många if förhållanden. De extra trådarna utför några onödiga tillägg, men vi bryr oss inte längre om innehållet de producerar. Eftersom varp körs i SIMD-läge sparar vi inte faktiskt tid genom att låta de trådarna göra ingenting. Å andra sidan tar det relativt lång tid att utvärdera förhållandena, eftersom dessa if uttalanden är så små är stora. Det ursprungliga if uttalet kan också tas bort om shArr[32..47] är stoppad med 0.
Minskningen av varpnivån kan också användas för att öka blocknivåreduktionen:
__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];
}
}
Argumentet &r[idx & ~(warpSize-1)] är i princip r + warpIdx*32 . Detta delar effektivt r arrayen i bitar av 32 element, och varje bit är tilldelad till separat varp.
Parallell minskning med en varp för icke-kommutativ operatör
Ibland måste reduktionen utföras i mycket liten skala, som en del av en större CUDA-kärna. Anta till exempel att inmatningsdata har exakt 32 element - antalet trådar i en varp. I ett sådant scenario kan ett enda varp tilldelas för att utföra reduktionen. Med tanke på att varp körs i en perfekt synk kan många __syncthreads() instruktioner tas bort - jämfört med en blocknivåreduktion.
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 är företrädesvis en matris i delat minne. Värdet ska vara detsamma för alla trådar i varpen. Om sumCommSingleWarp kallas av flera varp ska shArr vara olika mellan varp (samma inom varje varp).
Argumentet shArr markeras som volatile att säkerställa att operationer på arrayen faktiskt utförs där det anges. Annars kan den repetitiva tilldelningen till shArr[idx] optimeras som en tilldelning till ett register, där endast slutlig tilldelning är en verklig butik för shArr . När det händer är de omedelbara tilldelningarna inte synliga för andra trådar, vilket ger felaktiga resultat. Observera att du kan passera en normal icke-flyktig matris som ett argument för flyktig, samma som när du passerar non-const som en const-parameter.
Om man inte bryr sig om det slutliga innehållet i shArr[1..31] och kan shArr[32..47] med nollor, kan man förenkla koden ovan:
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];
}
I denna installation tog vi bort alla if förhållanden, som utgör ungefär hälften av instruktionerna. De extra trådarna utför några onödiga tillägg och lagrar resultatet i celler av shArr som i slutändan inte har någon inverkan på det slutliga resultatet. Eftersom varp körs i SIMD-läge sparar vi inte faktiskt tid genom att låta de trådarna göra ingenting.
Parallell reduktion med en varp med endast register
Normalt utförs reduktion på global eller delad matris. Men när reduktionen utförs i mycket liten skala, som en del av en större CUDA-kärna, kan den utföras med ett enda varp. När det händer, på Keppler eller högre arkitekturer (CC> = 3.0), är det möjligt att använda warp-shuffle-funktioner för att undvika att använda delat minne alls.
Anta till exempel att varje tråd i ett varp har ett enda inmatningsdatavärde. Alla trådar tillsammans har 32 element, som vi behöver sammanfatta (eller utföra annan associerande operation)
__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);
}
Denna version fungerar för både kommutativa och icke-kommutativa operatörer.