programacion - hola mundo en cuda

Encuentra el máximo de la matriz en CUDA

Question

Empecé en CUDA. Ahora tengo una pregunta. Tengo una matriz N * N, y una escala de ventana es 8x8. Quiero subdividir esta matriz en múltiples submatrices y encontrar el valor máximo de esto. Por ejemplo, si tengo una matriz 64 * 64, entonces tendré 8 matrices pequeñas con una escala 8 * 8 y descubriré 8 valores máximos. Finalmente guardo todos los valores máximos en una nueva matriz, pero su orden siempre cambia. Quiero encontrar una solución para mantenerlos en orden

__global__ void calculate_emax_kernel(float emap[],float emax[], int img_height, int img_width,int windows_size)
{
int x_index = blockIdx.x*blockDim.x+threadIdx.x;
int y_index = blockIdx.y*blockDim.y+threadIdx.y;
int num_row_block = img_height/windows_size;
int num_col_block = img_width/windows_size;
__shared__ float window_elements[256];
__shared__ int counter;
__shared__ int emax_count;
if (threadIdx.x == 0) emax_count = 0;
__syncthreads();
int index;
int emax_idx = 0;
if(y_index >= img_height|| x_index >= img_width) return;
for(int i = 0; i < num_row_block; i++)
{
for(int j = 0; j < num_col_block; j++)
{
counter = 0;
if(y_index >= i*windows_size && y_index < (i+1)*windows_size
&& x_index >= j*windows_size && x_index < (j+1)*windows_size)
{
int idx = y_index*img_height + x_index;
index = atomicAdd(&counter, 1);
window_elements[index] = emap[idx];
__syncthreads();
// reduction
unsigned int k = (windows_size*windows_size)/2;
while(k != 0)
{
if(index < k)
{
window_elements[index] = fmaxf(window_elements[index], window_elements[index+k]);
}
k /= 2;
}
if(index == 0)
{
emax[i*num_row_block+j] = window_elements[index];
}
}
__syncthreads();
}
__syncthreads();
}
__syncthreads();
}

Esta es mi configuración

void construct_emax(float *input,float *output, int img_height, int img_width)
{
int windows_size = 4;
float * d_input, * d_output;
cudaMalloc(&d_input, img_width*img_height*sizeof(float));
cudaMalloc(&d_output, img_width*img_height*sizeof(float));
cudaMemcpy(d_input, input, img_width*img_height*sizeof(float), cudaMemcpyHostToDevice);
dim3 blocksize(16,16);
dim3 gridsize;
gridsize.x=(img_width+blocksize.x-1)/blocksize.x;
gridsize.y=(img_height+blocksize.y-1)/blocksize.y;
calculate_emax_kernel<<<gridsize,blocksize>>>(d_input,d_output,img_height,img_width,windows_size);
}

Answer 1

Con CUDA, la reducción paralela es difícil; la reducción paralela segmentada es más complicada. Ahora lo está haciendo en 2-D, y su segmento / ventana es más pequeño que el bloque de subprocesos.

Para un tamaño de ventana grande, no creo que sea un problema. Podría usar un bloque de hilos para reducir una ventana. Por ejemplo, si tiene una ventana de 16x16, simplemente podría usar un bloque de hilos de 16x16. Si tiene un tamaño de ventana aún mayor, por ejemplo 64x64, aún podría usar un bloque de hilos de 16x16. Primero reduzca la ventana de 64x64 a 16x16 elementos durante la carga de datos, luego reduzca a 1 escalar dentro del bloque de hilos.

Para el tamaño de ventana más pequeño que el tamaño del bloque, tendrá que reducir varias ventanas por bloque de subprocesos para un mayor rendimiento. Podría usar su configuración actual de bloque / red, donde cada bloque de 256 hilos (16x16) es responsable de 16 ventanas de 4x4. Pero esto no será óptimo porque cada envoltura de 32 hilos está organizada en dos partes (2x16). Esto no es bueno para el acceso a la memoria global fusionada , y es difícil mapear una urdimbre de 2x16 a una o más ventanas de 4x4 para una reducción paralela eficiente.

Alternativamente, le sugiero que utilice el bloque de hilos 1-D con 256 hilos. Cada hilo m reduce una ventana de m . Entonces podría usar una cuadrícula 2-D para cubrir toda la imagen.

const int m = window_size;
dim3 blocksize(256);
dim3 gridsize((img_width+255)/256, (img_height+m-1)/m);

En la función kernel, podrías

1. reducir cada ventana de m x m a un vector de 1x m durante la carga de datos globales;
2. use el método de reducción de árboles para reducir el vector 1x m a un escalar.

El siguiente código es una demostración conceptual que funciona cuando m es una potencia de 2 y m <= 32 . Podrías modificarlo más para una arbitraria m y una mejor verificación de límites.

#include <assert.h>
#include <cuda.h>
#include <thrust/device_vector.h>
__global__ void calculate_emax_kernel(const float* input, float* output,
int height, int width, int win_size,
int out_width) {
const int tid = threadIdx.x;
const int i = blockIdx.y * win_size;
const int j = blockIdx.x * 256 + tid;
const int win_id = j % win_size;
__shared__ float smax[256];
float tmax = -1e20;
if (j < width) {
for (int tile = 0; tile < win_size; tile++) {
if (i + tile < height) {
tmax = max(tmax, input[(i + tile) * width + j]);
}
}
}
smax[tid] = tmax;
for (int shift = win_size / 2; shift > 0; shift /= 2) {
if (win_id < shift) {
smax[tid] = max(smax[tid], smax[tid + shift]);
}
}
if (win_id == 0 && j < width) {
output[blockIdx.y * out_width + (j / win_size)] = smax[tid];
}
}
int main() {
const int height = 1024;
const int width = 1024;
const int m = 4;
thrust::device_vector<float> in(height * width);
thrust::device_vector<float> out(
((height + m - 1) / m) * ((width + m - 1) / m));
dim3 blocksize(256);
dim3 gridsize((width + 255) / 256, (height + m - 1) / m);
assert(m == 2 || m == 4 || m == 8 || m == 16 || m == 32);
calculate_emax_kernel<<<gridsize, blocksize>>>(
thrust::raw_pointer_cast(in.data()),
thrust::raw_pointer_cast(out.data()),
height, width, m, (width + m - 1) / m);
return 0;
}