cuda
Parallelle reductie (bijv. Hoe een array op te tellen)

algorithm C Language C# Language C++ Java Language GNU/Linux MATLAB Language opencv Python Language tensorflow
Zoeken…


Opmerkingen

Parallelle reductie-algoritme verwijst meestal naar een algoritme dat een reeks elementen combineert en een enkel resultaat oplevert. Typische problemen die in deze categorie vallen zijn:

  • alle elementen in een array samenvatten
  • het vinden van een maximum in een array

Over het algemeen kan de parallelle reductie worden toegepast voor elke binaire associatieve operator , dwz (A*B)*C = A*(B*C) . Met een dergelijke operator * groepeert het parallelle reductie-algoritme de arrayargumenten in paren. Elk paar wordt parallel met anderen berekend, waardoor de totale reeksgrootte in één stap wordt gehalveerd. Het proces wordt herhaald totdat er slechts één element bestaat.

Als de operator commutatief is (dwz A*B = B*A ) en niet associatief is, kan het algoritme in een ander patroon worden gekoppeld. Vanuit theoretisch oogpunt maakt het geen verschil, maar in de praktijk geeft het een beter geheugentoegangspatroon:

Niet alle associatieve operatoren zijn commutatief - neem bijvoorbeeld matrixvermenigvuldiging.

Parallelle reductie met één blok voor commutatieve operator

De eenvoudigste aanpak voor parallelle reductie in CUDA is om een enkel blok toe te wijzen om de taak uit te voeren: 

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);

Dit is het meest haalbaar wanneer de gegevensgrootte niet erg groot is (ongeveer enkele duizenden elementen). Dit gebeurt meestal wanneer de reductie een onderdeel is van een groter CUDA-programma. Als de invoer vanaf het begin overeenkomt met blockSize , kan de eerste for lus volledig worden verwijderd.

Merk op dat in de eerste stap, wanneer er meer elementen zijn dan threads, we dingen volledig onafhankelijk optellen. Alleen wanneer het probleem is beperkt tot blockSize , wordt de daadwerkelijke parallelle reductie geactiveerd. Dezelfde code kan worden toegepast op elke andere commutatieve, associatieve operator, zoals vermenigvuldiging, minimum, maximum, enz.

Merk op dat het algoritme sneller kan worden gemaakt, bijvoorbeeld door een parallelle reductie op warp-niveau te gebruiken.

Parallelle reductie met één blok voor niet-commutatieve operator

Het uitvoeren van parallelle reductie voor een niet-commutatieve operator is iets ingewikkelder in vergelijking met de commutatieve versie. In het voorbeeld gebruiken we nog steeds een toevoeging over gehele getallen omwille van de eenvoud. Het zou bijvoorbeeld kunnen worden vervangen door matrixvermenigvuldiging die echt niet-commutatief is. Merk op dat 0 hierbij moet worden vervangen door een neutraal element van de vermenigvuldiging - dat wil zeggen een identiteitsmatrix. 

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);

In de eerste while-lus wordt uitgevoerd zolang er meer invoerelementen zijn dan threads. Bij elke iteratie wordt een enkele reductie uitgevoerd en wordt het resultaat gecomprimeerd in de eerste helft van de shArr reeks. De tweede helft wordt dan gevuld met nieuwe gegevens.

Zodra alle gegevens zijn geladen vanuit gArr , wordt de tweede lus uitgevoerd. Nu comprimeren we het resultaat niet langer (wat een extra __syncthreads() kost). In elke stap heeft de thread n toegang tot het 2*n -de actieve element en telt het op met 2*n+1 -de element:

Er zijn veel manieren om dit eenvoudige voorbeeld verder te optimaliseren, bijvoorbeeld door het verlagen van het warp-niveau en door conflicten met gedeelde geheugenbanken te verwijderen.

Multi-block parallelle reductie voor commutatieve operator

Multi-block benadering van parallelle reductie in CUDA vormt een extra uitdaging, vergeleken met single-block benadering, omdat blokken beperkt zijn in communicatie. Het idee is om elk blok een deel van de invoerarray te laten berekenen en vervolgens een laatste blok te hebben om alle gedeeltelijke resultaten samen te voegen. Om dat te doen kan men twee kernels lanceren, impliciet een roosterbreed synchronisatiepunt creëren. 

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;
}

Men wil idealiter voldoende blokken lanceren om alle multiprocessors op de GPU bij volledige bezetting te verzadigen. Het overschrijden van dit aantal - met name het starten van zoveel threads als er elementen in de array zijn - is contraproductief. Hierdoor verhoogt het de rekenkracht niet meer, maar voorkomt het gebruik van de zeer efficiënte eerste lus.

Het is ook mogelijk om hetzelfde resultaat te verkrijgen met behulp van een enkele kernel, met behulp van 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;
}

Merk op dat de kernel sneller gemaakt kan worden, bijvoorbeeld door een parallelle reductie op warp-niveau te gebruiken.

Multi-block parallelle reductie voor niet-commutatieve operator

Multi-block benadering van parallelle reductie lijkt erg op de single-block benadering. De globale invoerarray moet worden opgesplitst in secties, elk gereduceerd met een enkel blok. Wanneer een gedeeltelijk resultaat uit elk blok wordt verkregen, vermindert een laatste blok deze om het uiteindelijke resultaat te verkrijgen.

  • sumNoncommSingleBlock wordt meer in detail uitgelegd in het voorbeeld van reductie met één blok.
  • lastBlock accepteert alleen het laatste blok dat het bereikt. Als je dit wilt vermijden, kun je de kernel in twee afzonderlijke aanroepen splitsen. 
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);    
}

Men wil idealiter voldoende blokken lanceren om alle multiprocessors op de GPU bij volledige bezetting te verzadigen. Het overschrijden van dit aantal - met name het starten van zoveel threads als er elementen in de array zijn - is contraproductief. Hierdoor verhoogt het de rekenkracht niet meer, maar voorkomt het gebruik van de zeer efficiënte eerste lus.

Parallelle reductie met enkele kromming voor commutatieve operator

Soms moet de reductie op zeer kleine schaal worden uitgevoerd, als onderdeel van een grotere CUDA-kernel. Stel bijvoorbeeld dat de invoergegevens exact 32 elementen hebben - het aantal threads in een warp. In een dergelijk scenario kan een enkele kromming worden toegewezen om de reductie uit te voeren. Omdat warp perfect synchroon wordt uitgevoerd, kunnen veel __syncthreads() instructies worden verwijderd - in vergelijking met een reductie op blokniveau. 

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 is bij voorkeur een array in gedeeld geheugen. De waarde moet hetzelfde zijn voor alle threads in de ketting. Als sumCommSingleWarp door meerdere warps wordt shArr moet shArr tussen de warps verschillen (hetzelfde binnen elke warp).

Het argument shArr is gemarkeerd als volatile om ervoor te zorgen dat bewerkingen op de array daadwerkelijk worden uitgevoerd waar aangegeven. Anders kan de repetitieve toewijzing aan shArr[idx] worden geoptimaliseerd als een toewijzing aan een register, waarbij alleen de definitieve toewijzing een daadwerkelijke winkel voor shArr . Wanneer dat gebeurt, zijn de onmiddellijke toewijzingen niet zichtbaar voor andere threads, wat onjuiste resultaten oplevert. Merk op dat u een normale niet-vluchtige array kunt doorgeven als een argument van vluchtige, hetzelfde als wanneer u niet-const doorgeeft als const parameter.

Als men na de reductie niet shArr[1..31] de inhoud van shArr[1..31] , kan de code nog verder worden vereenvoudigd: 

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 deze opstelling hebben we veel if voorwaarden verwijderd. De extra threads voeren een aantal onnodige toevoegingen uit, maar we geven niet langer om de inhoud die ze produceren. Omdat warps worden uitgevoerd in de SIMD-modus, besparen we niet echt op tijd door die threads niets te laten doen. Aan de andere kant, het evalueren van de voorwaarden neemt relatief grote hoeveelheid tijd, omdat het lichaam van deze if uitspraken zijn zo klein. De initiële if instructie kan ook worden verwijderd als shArr[32..47] wordt opgevuld met 0.

De verlaging op kettingniveau kan ook worden gebruikt om de verlaging op blokniveau te stimuleren: 

__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];
    }
}

Het argument &r[idx & ~(warpSize-1)] is in principe r + warpIdx*32 . Dit splitst de r reeks effectief in brokken van 32 elementen, en elke brok wordt toegewezen aan afzonderlijke ketting.

Parallelle reductie met enkele kromming voor niet-commutatieve operator

Soms moet de reductie op zeer kleine schaal worden uitgevoerd, als onderdeel van een grotere CUDA-kernel. Stel bijvoorbeeld dat de invoergegevens exact 32 elementen hebben - het aantal threads in een warp. In een dergelijk scenario kan een enkele kromming worden toegewezen om de reductie uit te voeren. Omdat warp perfect synchroon wordt uitgevoerd, kunnen veel __syncthreads() instructies worden verwijderd - in vergelijking met een reductie op blokniveau. 

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 is bij voorkeur een array in gedeeld geheugen. De waarde moet hetzelfde zijn voor alle threads in de ketting. Als sumCommSingleWarp door meerdere warps wordt shArr moet shArr tussen de warps verschillen (hetzelfde binnen elke warp).

Het argument shArr is gemarkeerd als volatile om ervoor te zorgen dat bewerkingen op de array daadwerkelijk worden uitgevoerd waar aangegeven. Anders kan de repetitieve toewijzing aan shArr[idx] worden geoptimaliseerd als een toewijzing aan een register, waarbij alleen de definitieve toewijzing een daadwerkelijke winkel voor shArr . Wanneer dat gebeurt, zijn de onmiddellijke toewijzingen niet zichtbaar voor andere threads, wat onjuiste resultaten oplevert. Merk op dat u een normale niet-vluchtige array kunt doorgeven als een argument van vluchtige, hetzelfde als wanneer u niet-const doorgeeft als const parameter.

Als men niet shArr[1..31] de uiteindelijke inhoud van shArr[1..31] en shArr[32..47] met nullen kunt shArr[32..47] , kan men de bovenstaande code vereenvoudigen: 

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 deze opstelling hebben we alle if voorwaarden verwijderd, die ongeveer de helft van de instructies vormen. De extra threads voeren een aantal onnodige toevoegingen uit, waarbij het resultaat wordt opgeslagen in cellen van shArr die uiteindelijk geen invloed hebben op het eindresultaat. Omdat warps worden uitgevoerd in de SIMD-modus, besparen we niet echt op tijd door die threads niets te laten doen.

Parallelle reductie met enkele kromming, alleen met behulp van registers

Gewoonlijk wordt reductie uitgevoerd op globale of gedeelde array. Wanneer de reductie echter op zeer kleine schaal wordt uitgevoerd, als onderdeel van een grotere CUDA-kernel, kan deze worden uitgevoerd met een enkele warp. Wanneer dat gebeurt, op Keppler of hogere architecturen (CC> = 3.0), is het mogelijk om warp-shuffle-functies te gebruiken om helemaal geen gedeeld geheugen te gebruiken.

Stel bijvoorbeeld dat elke thread in een warp een enkele ingangsgegevenswaarde bevat. Alle threads samen hebben 32 elementen, die we moeten samenvatten (of een andere associatieve bewerking uitvoeren) 

__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);
}

Deze versie werkt voor zowel commutatieve als niet-commutatieve operatoren.



Modified text is an extract of the original Stack Overflow Documentation
Licentie onder CC BY-SA 3.0
Niet aangesloten bij Stack Overflow