cuda
Reducción paralela (por ejemplo, cómo sumar una matriz)
Buscar..
Observaciones
El algoritmo de reducción paralela generalmente se refiere a un algoritmo que combina una matriz de elementos, produciendo un único resultado. Los problemas típicos que entran en esta categoría son:
- resumiendo todos los elementos en una matriz
- encontrar un máximo en una matriz
En general, la reducción paralela se puede aplicar a cualquier operador asociativo binario, es decir (A*B)*C = A*(B*C) . Con tal operador *, el algoritmo de reducción paralelo agrupa repetidamente los argumentos de la matriz en pares. Cada par se calcula en paralelo con otros, reduciendo a la mitad el tamaño del arreglo en un solo paso. El proceso se repite hasta que solo existe un elemento.
Si el operador es conmutativo (es decir, A*B = B*A ) además de ser asociativo, el algoritmo puede emparejarse en un patrón diferente. Desde el punto de vista teórico, no hace ninguna diferencia, pero en la práctica proporciona un mejor patrón de acceso a la memoria:
No todos los operadores asociativos son conmutativos, como la multiplicación de matrices, por ejemplo.
Reducción en paralelo de bloque único para operador conmutativo
El enfoque más simple para la reducción paralela en CUDA es asignar un solo bloque para realizar la tarea:
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);
Esto es más factible cuando el tamaño de los datos no es muy grande (alrededor de unos pocos elementos). Esto suele ocurrir cuando la reducción es parte de un programa CUDA más grande. Si la entrada coincide con blockSize desde el principio, el primer bucle for se puede eliminar por completo.
Tenga en cuenta que en el primer paso, cuando hay más elementos que hilos, agregamos cosas de forma completamente independiente. Solo cuando el problema se reduce a blockSize , se activa la reducción paralela real. El mismo código se puede aplicar a cualquier otro operador asociativo, como conmutación, mínimo, máximo, etc.
Tenga en cuenta que el algoritmo se puede hacer más rápido, por ejemplo, mediante el uso de una reducción paralela de nivel de deformación.
Reducción en paralelo de bloque único para operador no conmutativo
Hacer una reducción paralela para un operador no conmutativo es un poco más complicado, en comparación con la versión conmutativa. En el ejemplo, todavía utilizamos una suma sobre enteros para simplificar. Podría reemplazarse, por ejemplo, con la multiplicación de matrices que realmente no es conmutativa. Tenga en cuenta que, al hacerlo, 0 debe reemplazarse por un elemento neutral de la multiplicación, es decir, una matriz de identidad.
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);
En el primer ciclo, el bucle se ejecuta siempre que haya más elementos de entrada que hilos. En cada iteración, se realiza una reducción única y el resultado se comprime en la primera mitad de la matriz shArr . La segunda mitad se llena con nuevos datos.
Una vez que todos los datos se cargan desde gArr , se ejecuta el segundo bucle. Ahora, ya no comprimimos el resultado (lo que cuesta __syncthreads() adicionales). En cada paso, el subproceso n accede al elemento activo 2*n -th y lo suma con el elemento 2*n+1 -th:
Hay muchas formas de optimizar aún más este simple ejemplo, por ejemplo, a través de la reducción del nivel de distorsión y eliminando los conflictos del banco de memoria compartida.
Reducción paralela multibloque para operador conmutativo.
El enfoque de bloques múltiples para la reducción paralela en CUDA plantea un desafío adicional, en comparación con el enfoque de un solo bloque, porque los bloques están limitados en la comunicación. La idea es dejar que cada bloque calcule una parte de la matriz de entrada y luego tener un bloque final para combinar todos los resultados parciales. Para hacer eso, uno puede lanzar dos núcleos, creando implícitamente un punto de sincronización en toda la red.
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;
}
Lo ideal sería lanzar bloques suficientes para saturar todos los multiprocesadores en la GPU en plena ocupación. Superar este número, en particular, lanzar tantos subprocesos como elementos en la matriz, es contraproducente. Al hacerlo, ya no aumenta la potencia de cálculo sin procesar, pero evita el uso del primer bucle muy eficiente.
También es posible obtener el mismo resultado utilizando un solo kernel, con la ayuda de la última guardia de bloque :
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;
}
Tenga en cuenta que el kernel se puede hacer más rápido, por ejemplo, utilizando una reducción paralela de nivel de deformación.
Reducción paralela multibloque para operador no conmutativo
El enfoque de bloques múltiples para la reducción paralela es muy similar al enfoque de un solo bloque. La matriz de entrada global debe dividirse en secciones, cada una reducida en un solo bloque. Cuando se obtiene un resultado parcial de cada bloque, un bloque final los reduce para obtener el resultado final.
-
sumNoncommSingleBlockse explica con más detalle en el ejemplo de reducción de un solo bloque. -
lastBlockacepta solo el último bloque que llega a él. Si quiere evitar esto, puede dividir el kernel en dos llamadas separadas.
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);
}
Lo ideal sería lanzar bloques suficientes para saturar todos los multiprocesadores en la GPU en plena ocupación. Superar este número, en particular, lanzar tantos subprocesos como elementos en la matriz, es contraproducente. Al hacerlo, ya no aumenta la potencia de cálculo sin procesar, pero evita el uso del primer bucle muy eficiente.
Reducción paralela de una sola urdimbre para operador conmutativo
A veces, la reducción debe realizarse en una escala muy pequeña, como parte de un núcleo CUDA más grande. Supongamos, por ejemplo, que los datos de entrada tienen exactamente 32 elementos: el número de subprocesos en una deformación. En tal escenario, se puede asignar una sola urdimbre para realizar la reducción. Dado que la deformación se ejecuta en una sincronización perfecta, muchas instrucciones __syncthreads() se pueden eliminar, en comparación con una reducción a nivel de bloque.
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 es preferiblemente una matriz en memoria compartida. El valor debe ser el mismo para todos los subprocesos en la deformación. Si sumCommSingleWarp es llamado por múltiples hilos de urdimbre, shArr debe ser diferente entre urdimbres (las mismas dentro de cada urdimbre).
El argumento shArr se marca como volatile para garantizar que las operaciones en la matriz se realicen donde se indica. De lo contrario, la asignación shArr[idx] a shArr[idx] se puede optimizar como una asignación a un registro, y solo la asignación final es una tienda real para shArr . Cuando eso sucede, las asignaciones inmediatas no son visibles para otros subprocesos, lo que produce resultados incorrectos. Tenga en cuenta que puede pasar una matriz normal no volátil como un argumento de volatile, igual que cuando pasa no constante como un parámetro const.
Si a uno no le importa el contenido de shArr[1..31] después de la reducción, puede simplificar el código aún más:
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];
}
En esta configuración eliminamos muchas if condiciones. Los hilos adicionales realizan algunas adiciones innecesarias, pero ya no nos preocupamos por los contenidos que producen. Dado que las deformaciones se ejecutan en modo SIMD, en realidad no ahorramos a tiempo porque esos subprocesos no hacen nada. Por otro lado, la evaluación de las condiciones lleva relativamente mucho tiempo, ya que el cuerpo de estas if declaraciones son tan pequeñas. La sentencia if inicial también se puede eliminar si shArr[32..47] se rellena con 0.
La reducción de nivel de deformación se puede utilizar para aumentar la reducción de nivel de bloque también:
__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];
}
}
El argumento &r[idx & ~(warpSize-1)] es básicamente r + warpIdx*32 . Esto divide efectivamente la matriz r en trozos de 32 elementos, y cada trozo se asigna a warp separado.
Reducción paralela de una sola urdimbre para un operador no conmutativo.
A veces, la reducción debe realizarse en una escala muy pequeña, como parte de un núcleo CUDA más grande. Supongamos, por ejemplo, que los datos de entrada tienen exactamente 32 elementos: el número de subprocesos en una deformación. En tal escenario, se puede asignar una sola urdimbre para realizar la reducción. Dado que la deformación se ejecuta en una sincronización perfecta, muchas instrucciones __syncthreads() se pueden eliminar, en comparación con una reducción a nivel de bloque.
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 es preferiblemente una matriz en memoria compartida. El valor debe ser el mismo para todos los subprocesos en la deformación. Si sumCommSingleWarp es llamado por múltiples hilos de urdimbre, shArr debe ser diferente entre urdimbres (las mismas dentro de cada urdimbre).
El argumento shArr se marca como volatile para garantizar que las operaciones en la matriz se realicen donde se indica. De lo contrario, la asignación shArr[idx] a shArr[idx] se puede optimizar como una asignación a un registro, y solo la asignación final es una tienda real para shArr . Cuando eso sucede, las asignaciones inmediatas no son visibles para otros subprocesos, lo que produce resultados incorrectos. Tenga en cuenta que puede pasar una matriz normal no volátil como un argumento de volatile, igual que cuando pasa no constante como un parámetro const.
Si a uno no le importa el contenido final de shArr[1..31] y puede rellenar shArr[32..47] con ceros, puede simplificar el código anterior:
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];
}
En esta configuración hemos eliminado todo if las condiciones, que constituyen aproximadamente la mitad de las instrucciones. Los subprocesos adicionales realizan algunas adiciones innecesarias, almacenando el resultado en celdas de shArr que, en última instancia, no tienen impacto en el resultado final. Dado que las deformaciones se ejecutan en modo SIMD, en realidad no ahorramos a tiempo porque esos subprocesos no hacen nada.
Reducción paralela de una sola urdimbre usando solo registros
Normalmente, la reducción se realiza en una matriz global o compartida. Sin embargo, cuando la reducción se realiza en una escala muy pequeña, como parte de un núcleo CUDA más grande, se puede realizar con una sola deformación. Cuando eso sucede, en las arquitecturas Keppler o superiores (CC> = 3.0), es posible usar las funciones de warp-shuffle para evitar el uso de la memoria compartida.
Supongamos, por ejemplo, que cada subproceso en una deformación contiene un solo valor de datos de entrada. Todos los hilos juntos tienen 32 elementos, que necesitamos resumir (o realizar otra operación asociativa)
__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);
}
Esta versión funciona para operadores conmutativos y no conmutativos.