cuda
Parallele Reduktion (z. B. wie ein Array summiert wird)
Suche…
Bemerkungen
Paralleler Reduktionsalgorithmus bezieht sich typischerweise auf einen Algorithmus, der ein Array von Elementen kombiniert, um ein einzelnes Ergebnis zu erzeugen. Typische Probleme, die in diese Kategorie fallen, sind:
- Aufsummieren aller Elemente in einem Array
- ein Maximum in einem Array finden
Im Allgemeinen kann die parallele Reduktion für jeden binären assoziativen Operator angewendet werden , dh (A*B)*C = A*(B*C) . Mit einem solchen Operator * gruppiert der Parallelreduktionsalgorithmus die Array-Argumente wiederholt in Paaren. Jedes Paar wird parallel zu anderen berechnet, wobei die Gesamtgröße des Arrays in einem Schritt halbiert wird. Der Vorgang wird wiederholt, bis nur ein einziges Element vorhanden ist.
Wenn der Operator kommutativ ist (dh A*B = B*A ) und nicht assoziativ ist, kann der Algorithmus in einem anderen Muster paaren. Vom theoretischen Standpunkt aus macht es keinen Unterschied, aber in der Praxis ergibt sich ein besseres Speicherzugriffsmuster:
Nicht alle assoziativen Operatoren sind kommutativ - nehmen Sie zum Beispiel die Matrixmultiplikation.
Parallele Reduktion um einen Block für kommutative Operatoren
Der einfachste Ansatz für die parallele Reduzierung von CUDA besteht darin, einen einzelnen Block zur Ausführung der Aufgabe zuzuweisen:
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);
Dies ist am besten machbar, wenn die Datengröße nicht sehr groß ist (etwa einige tausend Elemente). Dies geschieht normalerweise, wenn die Reduktion Teil eines größeren CUDA-Programms ist. Wenn die Eingabe von Anfang an mit blockSize , kann die erste for Schleife vollständig entfernt werden.
Beachten Sie, dass wir im ersten Schritt, wenn es mehr Elemente als Threads gibt, die Dinge vollständig unabhängig voneinander addieren. Nur wenn das Problem auf blockSize reduziert blockSize , wird die tatsächliche parallele Reduktion blockSize . Derselbe Code kann auf jeden anderen kommutativen, assoziativen Operator angewendet werden, wie Multiplikation, Minimum, Maximum usw.
Beachten Sie, dass der Algorithmus schneller gemacht werden kann, z. B. durch eine parallele Reduktion auf Warp-Ebene.
Parallele Reduktion um einen Block für nicht kommutativen Operator
Die parallele Reduzierung für einen nicht kommutativen Operator ist im Vergleich zur kommutativen Version etwas komplizierter. Im Beispiel verwenden wir der Einfachheit halber immer noch eine Addition über Ganzzahlen. Sie könnte beispielsweise durch eine Matrixmultiplikation ersetzt werden, die eigentlich nicht kommutativ ist. Beachten Sie dabei, dass 0 durch ein neutrales Element der Multiplikation ersetzt werden sollte, dh eine Identitätsmatrix.
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 der ersten while-Schleife wird ausgeführt, solange mehr Eingabeelemente als Threads vorhanden sind. Bei jeder Iteration wird eine einzelne Reduktion durchgeführt und das Ergebnis in die erste Hälfte des shArr Arrays shArr . Die zweite Hälfte wird dann mit neuen Daten gefüllt.
Sobald alle Daten von gArr geladen gArr , wird die zweite Schleife ausgeführt. Jetzt komprimieren wir das Ergebnis nicht mehr (was __syncthreads() zusätzlich __syncthreads() ). In jedem Schritt greift der Thread n auf das 2*n te aktive Element zu und summiert es mit 2*n+1 ten Element auf:
Es gibt viele Möglichkeiten, dieses einfache Beispiel weiter zu optimieren, z. B. durch Reduzierung des Warp-Pegels und durch das Entfernen von Konflikten in der gemeinsamen Speicherbank.
Parallele Multi-Block-Reduktion für kommutativen Operator
Ein Ansatz mit mehreren Blöcken zur parallelen Reduktion von CUDA stellt im Vergleich zum Ansatz mit einem Block eine zusätzliche Herausforderung dar, da Blöcke in der Kommunikation begrenzt sind. Die Idee ist, dass jeder Block einen Teil des Eingabearrays berechnen lässt und dann einen letzten Block hat, um alle Teilergebnisse zusammenzuführen. Dazu kann man zwei Kernel starten und implizit einen gitterweiten Synchronisationspunkt erstellen.
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;
}
Idealerweise möchte man genügend Blöcke starten, um alle Multiprozessoren der GPU bei voller Belegung zu sättigen. Das Überschreiten dieser Anzahl, insbesondere das Starten von so vielen Threads, wie Elemente im Array vorhanden sind, ist kontraproduktiv. Dadurch wird die reine Rechenleistung nicht mehr erhöht, sondern die Verwendung der sehr effizienten ersten Schleife verhindert.
Es ist auch möglich, mit Hilfe des Last-Block-Guard das gleiche Ergebnis mit einem einzelnen Kernel zu erzielen:
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;
}
Beachten Sie, dass der Kernel schneller gemacht werden kann, beispielsweise durch eine parallele Reduktion auf Warp-Ebene.
Parallele Multi-Block-Reduktion für nichtkommutativen Operator
Der Multi-Block-Ansatz für die parallele Reduktion ist dem Single-Block-Ansatz sehr ähnlich. Das globale Eingabearray muss in Abschnitte unterteilt werden, die jeweils um einen einzelnen Block reduziert sind. Wenn ein Teilergebnis von jedem Block erhalten wird, reduziert ein letzter Block diese, um das Endergebnis zu erhalten.
-
sumNoncommSingleBlockwird im Beispiel zur Reduzierung von Einzelblöcken ausführlicher erklärt. -
lastBlockakzeptiert nur den letzten Block, der ihn erreicht. Wenn Sie dies vermeiden möchten, können Sie den Kernel in zwei separate Aufrufe aufteilen.
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);
}
Idealerweise möchte man genügend Blöcke starten, um alle Multiprozessoren der GPU bei voller Belegung zu sättigen. Das Überschreiten dieser Anzahl, insbesondere das Starten von so vielen Threads, wie Elemente im Array vorhanden sind, ist kontraproduktiv. Dadurch wird die reine Rechenleistung nicht mehr erhöht, sondern die Verwendung der sehr effizienten ersten Schleife verhindert.
Single-Warp-Parallelreduktion für kommutative Operatoren
Manchmal muss die Reduktion in sehr kleinem Umfang als Teil eines größeren CUDA-Kernels durchgeführt werden. Angenommen, die Eingabedaten haben genau 32 Elemente - die Anzahl der Threads in einem Warp. In einem solchen Szenario kann ein einzelner Warp zugewiesen werden, um die Reduzierung durchzuführen. Da Warp in einer perfekten Synchronisation ausgeführt wird, können viele __syncthreads() Anweisungen entfernt werden - verglichen mit einer Reduktion auf Blockebene.
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 ist vorzugsweise ein Array im Shared Memory. Der Wert sollte für alle Threads im Warp derselbe sein. Wenn sumCommSingleWarp von mehreren Warps shArr sollte shArr zwischen Warps unterschiedlich sein (in jedem Warp gleich).
Das Argument shArr wird als volatile markiert, um sicherzustellen, dass Operationen an dem Array tatsächlich ausgeführt werden, shArr angegeben. Andernfalls kann die wiederholte Zuordnung zu shArr[idx] als Zuordnung zu einem Register optimiert werden, wobei nur die endgültige Zuordnung ein tatsächlicher Speicher für shArr . In diesem Fall sind die unmittelbaren Zuweisungen für andere Threads nicht sichtbar, was zu falschen Ergebnissen führt. Beachten Sie, dass Sie ein normales nichtflüchtiges Array als Argument eines flüchtigen Arrays übergeben können, genauso wie wenn Sie non-const als const-Parameter übergeben.
Wenn sich der Inhalt von shArr[1..31] nach der Reduktion nicht interessiert, kann der Code noch weiter vereinfacht werden:
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 diesem Setup haben wir viele if Bedingungen entfernt. Die zusätzlichen Threads führen einige unnötige Ergänzungen aus, aber wir kümmern uns nicht mehr um die Inhalte, die sie produzieren. Da Warps im SIMD-Modus ausgeführt werden, sparen wir eigentlich keine Zeit, da diese Threads nichts unternehmen. Auf der anderen Seite dauert die Bewertung der Bedingungen relativ viel Zeit, da der Rumpf dieser if Anweisungen so klein ist. Die anfängliche if Anweisung kann auch entfernt werden, wenn shArr[32..47] mit 0 aufgefüllt wird.
Die Reduzierung des Warp-Pegels kann auch verwendet werden, um die Reduktion auf Blockebene zu beschleunigen:
__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];
}
}
Das Argument &r[idx & ~(warpSize-1)] lautet grundsätzlich r + warpIdx*32 . Dadurch wird das r Array effektiv in Blöcke mit 32 Elementen aufgeteilt, und jeder Block wird einem separaten Warp zugewiesen.
Single-Warp-Parallelreduktion für nichtkommutative Operatoren
Manchmal muss die Reduktion in sehr kleinem Umfang als Teil eines größeren CUDA-Kernels durchgeführt werden. Angenommen, die Eingabedaten haben genau 32 Elemente - die Anzahl der Threads in einem Warp. In einem solchen Szenario kann ein einzelner Warp zugewiesen werden, um die Reduzierung durchzuführen. Da Warp in einer perfekten Synchronisation ausgeführt wird, können viele __syncthreads() Anweisungen entfernt werden - verglichen mit einer Reduktion auf Blockebene.
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 ist vorzugsweise ein Array im Shared Memory. Der Wert sollte für alle Threads im Warp derselbe sein. Wenn sumCommSingleWarp von mehreren Warps shArr sollte shArr zwischen Warps unterschiedlich sein (in jedem Warp gleich).
Das Argument shArr wird als volatile markiert, um sicherzustellen, dass Operationen an dem Array tatsächlich ausgeführt werden, shArr angegeben. Andernfalls kann die wiederholte Zuordnung zu shArr[idx] als Zuordnung zu einem Register optimiert werden, wobei nur die endgültige Zuordnung ein tatsächlicher Speicher für shArr . In diesem Fall sind die unmittelbaren Zuweisungen für andere Threads nicht sichtbar, was zu falschen Ergebnissen führt. Beachten Sie, dass Sie ein normales nichtflüchtiges Array als Argument eines flüchtigen Arrays übergeben können, genauso wie wenn Sie non-const als const-Parameter übergeben.
Wenn man sich nicht für den endgültigen Inhalt von shArr[1..31] und shArr[32..47] mit Nullen shArr[32..47] kann, kann man den obigen Code vereinfachen:
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 diesem Setup haben wir alle if Bedingungen entfernt, die etwa die Hälfte der Anweisungen ausmachen. Die zusätzlichen Threads führen einige unnötige Hinzufügungen durch und speichern das Ergebnis in Zellen von shArr , die letztendlich keinen Einfluss auf das Endergebnis haben. Da Warps im SIMD-Modus ausgeführt werden, sparen wir eigentlich keine Zeit, da diese Threads nichts unternehmen.
Single-Warp-Parallel-Reduktion nur mit Registern
Normalerweise wird die Reduktion für globale oder gemeinsam genutzte Arrays durchgeführt. Wenn die Reduktion jedoch in sehr kleinem Umfang als Teil eines größeren CUDA-Kernels durchgeführt wird, kann sie mit einem einzelnen Warp durchgeführt werden. In diesem Fall ist es bei Keppler- oder höheren Architekturen (CC> = 3.0) möglich, Warp-Shuffle-Funktionen zu verwenden, um den gemeinsamen Speicher zu vermeiden.
Angenommen, jeder Thread in einem Warp enthält einen einzelnen Eingabedatenwert. Alle Threads zusammen haben 32 Elemente, die wir zusammenfassen müssen (oder andere assoziative Operationen ausführen).
__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);
}
Diese Version funktioniert sowohl für kommutative als auch für nicht kommutative Operatoren.