cuda
Redukcja równoległa (np. Jak zsumować tablicę)
Szukaj…
Uwagi
Algorytm równoległej redukcji zazwyczaj odnosi się do algorytmu, który łączy w sobie tablicę elementów, tworząc pojedynczy wynik. Typowe problemy należące do tej kategorii to:
- sumując wszystkie elementy w tablicy
- znajdowanie maksimum w tablicy
Zasadniczo redukcję równoległą można zastosować dla dowolnego binarnego operatora asocjacyjnego , tj. (A*B)*C = A*(B*C) . W przypadku takiego operatora * algorytm równoległej redukcji wielokrotnie grupuje argumenty tablic w parach. Każda para jest obliczana równolegle z innymi, zmniejszając o połowę ogólny rozmiar tablicy w jednym kroku. Proces powtarza się, aż istnieje tylko jeden element.
Jeśli operator jest przemienny (tj. A*B = B*A ) oprócz tego, że jest asocjatywny, algorytm może parować według innego wzorca. Z teoretycznego punktu widzenia nie ma to znaczenia, ale w praktyce daje lepszy wzór dostępu do pamięci:
Nie wszystkie operatory asocjacyjne są przemienne - weźmy na przykład mnożenie macierzy.
Redukcja równoległa w jednym bloku dla operatora przemiennego
Najprostszym podejściem do równoległej redukcji w CUDA jest przypisanie jednego bloku do wykonania zadania:
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);
Jest to najbardziej wykonalne, gdy rozmiar danych nie jest bardzo duży (około kilku tysięcy elementów). Zwykle dzieje się tak, gdy redukcja jest częścią większego programu CUDA. Jeśli dane wejściowe pasują do blockSize od samego początku, pierwszą pętlę for można całkowicie usunąć.
Zauważ, że w pierwszym kroku, gdy jest więcej elementów niż wątków, dodajemy wszystko całkowicie niezależnie. Dopiero gdy problem zostanie zredukowany do blockSize , blockSize się faktyczna równoległa redukcja. Ten sam kod można zastosować do dowolnego innego operatora przemiennego, asocjacyjnego, takiego jak mnożenie, minimum, maksimum itp.
Zauważ, że algorytm może być przyspieszony, na przykład poprzez zastosowanie równoległej redukcji równoległej na poziomie warp.
Redukcja równoległa w jednym bloku dla operatora nieprzemiennego
Redukcja równoległa dla operatora nieprzemiennego jest nieco bardziej zaangażowana w porównaniu z wersją przemienną. W tym przykładzie nadal używamy dodawania ponad liczbami całkowitymi dla uproszczenia. Można go na przykład zastąpić mnożeniem macierzy, co tak naprawdę nie jest przemienne. Zauważ, że robiąc to, 0 powinno być zastąpione neutralnym elementem mnożenia - tj. Macierzą tożsamości.
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);
W pierwszej chwili pętla wykonuje się tak długo, jak długo jest więcej elementów wejściowych niż wątków. W każdej iteracji wykonywana jest pojedyncza redukcja, a wynik jest kompresowany do pierwszej połowy tablicy shArr . Druga połowa jest następnie wypełniana nowymi danymi.
Po załadowaniu wszystkich danych z gArr , druga pętla jest wykonywana. Teraz nie kompresujemy już wyniku (co kosztuje dodatkowe __syncthreads() ). W każdym kroku wątek n uzyskuje dostęp do 2*n aktywnego elementu i dodaje go z 2*n+1 elementem:
Istnieje wiele sposobów dalszej optymalizacji tego prostego przykładu, np. Poprzez zmniejszenie poziomu wypaczenia i usunięcie konfliktów banków pamięci współdzielonej.
Wieloblokowa redukcja równoległa dla operatora przemiennego
Podejście z wieloma blokami do równoległej redukcji w CUDA stanowi dodatkowe wyzwanie w porównaniu z podejściem z jednym blokiem, ponieważ bloki mają ograniczoną komunikację. Chodzi o to, aby każdy blok mógł obliczyć część tablicy wejściowej, a następnie mieć jeden blok końcowy, aby scalić wszystkie wyniki częściowe. W tym celu można uruchomić dwa jądra, domyślnie tworząc punkt synchronizacji dla całej siatki.
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;
}
Najlepiej jest uruchomić wystarczającą liczbę bloków, aby nasycić wszystkie procesory wieloprocesorowe na GPU przy pełnym obciążeniu. Przekroczenie tej liczby - w szczególności uruchomienie tyle wątków, ile jest elementów w tablicy - przynosi efekt przeciwny do zamierzonego. Robiąc to, nie zwiększa już surowej mocy obliczeniowej, ale uniemożliwia użycie bardzo wydajnej pierwszej pętli.
Możliwe jest również uzyskanie tego samego wyniku za pomocą jednego jądra, przy pomocy zabezpieczenia ostatniego bloku :
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;
}
Zauważ, że jądro może być przyspieszone, na przykład poprzez zastosowanie równoległej redukcji równoległej na poziomie warp.
Wieloblokowa redukcja równoległa dla operatora nieprzemiennego
Podejście z wieloma blokami do równoległej redukcji jest bardzo podobne do podejścia z jednym blokiem. Globalna tablica wejściowa musi być podzielona na sekcje, każda zmniejszona o jeden blok. Po uzyskaniu częściowego wyniku z każdego bloku, jeden blok końcowy zmniejsza je, aby uzyskać wynik końcowy.
-
sumNoncommSingleBlockwyjaśniono bardziej szczegółowo w przykładzie redukcji pojedynczego bloku. -
lastBlockakceptuje tylko ostatni blok, który do niego dotarł. Jeśli chcesz tego uniknąć, możesz podzielić jądro na dwa osobne wywołania.
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);
}
Najlepiej jest uruchomić wystarczającą liczbę bloków, aby nasycić wszystkie procesory wieloprocesorowe na GPU przy pełnym obciążeniu. Przekroczenie tej liczby - w szczególności uruchomienie tyle wątków, ile jest elementów w tablicy - przynosi efekt przeciwny do zamierzonego. Robiąc to, nie zwiększa już surowej mocy obliczeniowej, ale uniemożliwia użycie bardzo wydajnej pierwszej pętli.
Redukcja równoległa pojedynczej osnowy dla operatora przemiennego
Czasami redukcję należy wykonać na bardzo małą skalę, jako część większego jądra CUDA. Załóżmy na przykład, że dane wejściowe mają dokładnie 32 elementy - liczbę wątków w warp. W takim scenariuszu można przypisać pojedynczą warp w celu przeprowadzenia redukcji. Biorąc pod uwagę, że wypaczanie wykonuje się w idealnej synchronizacji, wiele instrukcji __syncthreads() można usunąć - w porównaniu z redukcją na poziomie bloku.
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 jest najlepiej tablicą we wspólnej pamięci. Wartość powinna być taka sama dla wszystkich wątków w osnowie. Jeśli sumCommSingleWarp jest wywoływana przez wiele shArr , shArr powinien różnić się między wypaczeniami (tak samo w każdej wypaczeniu).
Argument shArr jest oznaczony jako volatile aby zapewnić, że operacje na tablicy są faktycznie wykonywane tam, gdzie jest to wskazane. W przeciwnym razie powtarzające się przypisanie do shArr[idx] może zostać zoptymalizowane jako przypisanie do rejestru, przy czym tylko ostateczne przypisanie jest faktycznym shArr dla shArr . Kiedy tak się dzieje, natychmiastowe przypisania nie są widoczne dla innych wątków, co daje nieprawidłowe wyniki. Zauważ, że możesz przekazać normalną nieulotną tablicę jako argument zmiennej niestabilnej, tak samo jak w przypadku podania zmiennej non-const jako parametru const.
Jeśli nie zależy nam na zawartości shArr[1..31] po redukcji, można jeszcze bardziej uprościć kod:
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];
}
W tej konfiguracji usunęliśmy wiele warunków if . Dodatkowe wątki wprowadzają niepotrzebne dodatki, ale nie dbamy już o ich treść. Ponieważ osnowy działają w trybie SIMD, tak naprawdę nie oszczędzamy na czasie, ponieważ nici nie robią nic. Z drugiej strony, ocena warunków zajmuje stosunkowo dużo czasu, ponieważ ich treść, if stwierdzenia są tak małe. Początkową instrukcję if można również usunąć, jeśli shArr[32..47] zostanie shArr[32..47] 0.
Redukcję poziomu osnowy można również wykorzystać do zwiększenia redukcji poziomu bloku:
__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];
}
}
Argument &r[idx & ~(warpSize-1)] to w zasadzie r + warpIdx*32 . To skutecznie dzieli tablicę r na fragmenty po 32 elementy, a każdy fragment jest przypisany do osobnego osnowy.
Redukcja równoległa pojedynczej osnowy dla nieprzemiennego operatora
Czasami redukcję należy wykonać na bardzo małą skalę, jako część większego jądra CUDA. Załóżmy na przykład, że dane wejściowe mają dokładnie 32 elementy - liczbę wątków w warp. W takim scenariuszu można przypisać pojedynczą warp w celu przeprowadzenia redukcji. Biorąc pod uwagę, że wypaczanie wykonuje się w idealnej synchronizacji, wiele instrukcji __syncthreads() można usunąć - w porównaniu z redukcją na poziomie bloku.
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 jest najlepiej tablicą we wspólnej pamięci. Wartość powinna być taka sama dla wszystkich wątków w osnowie. Jeśli sumCommSingleWarp jest wywoływana przez wiele shArr , shArr powinien różnić się między wypaczeniami (tak samo w każdej wypaczeniu).
Argument shArr jest oznaczony jako volatile aby zapewnić, że operacje na tablicy są faktycznie wykonywane tam, gdzie jest to wskazane. W przeciwnym razie powtarzające się przypisanie do shArr[idx] może zostać zoptymalizowane jako przypisanie do rejestru, przy czym tylko ostateczne przypisanie jest faktycznym shArr dla shArr . Kiedy tak się dzieje, natychmiastowe przypisania nie są widoczne dla innych wątków, co daje nieprawidłowe wyniki. Zauważ, że możesz przekazać normalną nieulotną tablicę jako argument zmiennej niestabilnej, tak samo jak w przypadku podania zmiennej non-const jako parametru const.
Jeśli ktoś nie dba o ostateczną zawartość shArr[1..31] i może shArr[32..47] zerami, można uprościć powyższy kod:
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];
}
W tej konfiguracji usunęliśmy wszystko if warunki, które stanowią około połowy instrukcji. Dodatkowe wątki dokonują niepotrzebnych dodatków, przechowując wynik w komórkach shArr które ostatecznie nie mają wpływu na końcowy wynik. Ponieważ osnowy działają w trybie SIMD, tak naprawdę nie oszczędzamy na czasie, ponieważ nici nie robią nic.
Redukcja równoległa pojedynczej osnowy przy użyciu tylko rejestrów
Zwykle redukcję przeprowadza się na tablicy globalnej lub współdzielonej. Jednak gdy redukcja jest wykonywana na bardzo małą skalę, jako część większego jądra CUDA, można ją wykonać za pomocą pojedynczej warp. Kiedy tak się dzieje, na architekturze Keppler lub wyższej (CC> = 3.0) można użyć funkcji tasowania warp, aby w ogóle uniknąć używania pamięci współdzielonej.
Załóżmy na przykład, że każdy wątek w osnowie zawiera jedną wartość danych wejściowych. Wszystkie wątki razem mają 32 elementy, które musimy podsumować (lub wykonać inną operację asocjacyjną)
__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);
}
Ta wersja działa zarówno dla operatorów przemiennych, jak i nieprzemiennych.