CUDA: memoria compartida en una matriz bidimensional 2D

Question

Tuve un problema CUDA simple para una asignación de clase, pero el profesor agregó una tarea opcional para implementar el mismo algoritmo utilizando memoria compartida. No pude terminarlo antes de la fecha límite (como en, la fecha de entrada fue hace una semana) pero todavía tengo curiosidad, así que ahora voy a preguntar en Internet;).

La tarea básica era implementar una versión bastardizada de una sucesiva relajación excesiva rojo-negro tanto secuencialmente como en CUDA, asegúrese de obtener el mismo resultado en ambos y luego compare la aceleración. Como dije, hacerlo con memoria compartida era un complemento opcional de + 10%.

Voy a publicar mi versión de trabajo y el pseudocódigo de lo que he intentado hacer ya que no tengo el código en mis manos en este momento, pero puedo actualizarlo más tarde con el código real si alguien necesita eso.

Antes de que nadie lo diga: Sí, sé que usar CUtil es poco convincente, pero hizo que la comparación y los temporizadores fueran más fáciles.

Trabajando versión de la memoria global:

#include <stdlib.h> 
#include <stdio.h> 
#include <cutil_inline.h> 

#define N 1024 

__global__ void kernel(int *d_A, int *d_B) { 
    unsigned int index_x = blockIdx.x * blockDim.x + threadIdx.x; 
    unsigned int index_y = blockIdx.y * blockDim.y + threadIdx.y; 

    // map the two 2D indices to a single linear, 1D index 
    unsigned int grid_width = gridDim.x * blockDim.x; 
    unsigned int index = index_y * grid_width + index_x; 

    // check for boundaries and write out the result 
    if((index_x > 0) && (index_y > 0) && (index_x < N-1) && (index_y < N-1)) 
     d_B[index] = (d_A[index-1]+d_A[index+1]+d_A[index+N]+d_A[index-N])/4; 

} 

main (int argc, char **argv) { 

    int A[N][N], B[N][N]; 
    int *d_A, *d_B; // These are the copies of A and B on the GPU 
    int *h_B; // This is a host copy of the output of B from the GPU 
    int i, j; 
    int num_bytes = N * N * sizeof(int); 

    // Input is randomly generated 
    for(i=0;i<N;i++) { 
     for(j=0;j<N;j++) { 
      A[i][j] = rand()/1795831; 
      //printf("%d\n",A[i][j]); 
     } 
    } 

    cudaEvent_t start_event0, stop_event0; 
    float elapsed_time0; 
    CUDA_SAFE_CALL(cudaEventCreate(&start_event0)); 
    CUDA_SAFE_CALL(cudaEventCreate(&stop_event0)); 
    cudaEventRecord(start_event0, 0); 
    // sequential implementation of main computation 
    for(i=1;i<N-1;i++) { 
     for(j=1;j<N-1;j++) { 
      B[i][j] = (A[i-1][j]+A[i+1][j]+A[i][j-1]+A[i][j+1])/4; 
     } 
    } 
    cudaEventRecord(stop_event0, 0); 
    cudaEventSynchronize(stop_event0); 
    CUDA_SAFE_CALL(cudaEventElapsedTime(&elapsed_time0,start_event0, stop_event0)); 



    h_B = (int *)malloc(num_bytes); 
    memset(h_B, 0, num_bytes); 
    //ALLOCATE MEMORY FOR GPU COPIES OF A AND B 
    cudaMalloc((void**)&d_A, num_bytes); 
    cudaMalloc((void**)&d_B, num_bytes); 
    cudaMemset(d_A, 0, num_bytes); 
    cudaMemset(d_B, 0, num_bytes); 

    //COPY A TO GPU 
    cudaMemcpy(d_A, A, num_bytes, cudaMemcpyHostToDevice); 

    // create CUDA event handles for timing purposes 
    cudaEvent_t start_event, stop_event; 
    float elapsed_time; 
    CUDA_SAFE_CALL(cudaEventCreate(&start_event)); 
    CUDA_SAFE_CALL(cudaEventCreate(&stop_event)); 
    cudaEventRecord(start_event, 0); 

// TODO: CREATE BLOCKS AND THREADS AND INVOKE GPU KERNEL 
    dim3 block_size(256,1,1); //values experimentally determined to be fastest 

    dim3 grid_size; 
    grid_size.x = N/block_size.x; 
    grid_size.y = N/block_size.y; 

    kernel<<<grid_size,block_size>>>(d_A,d_B); 

    cudaEventRecord(stop_event, 0); 
    cudaEventSynchronize(stop_event); 
    CUDA_SAFE_CALL(cudaEventElapsedTime(&elapsed_time,start_event, stop_event)); 

    //COPY B BACK FROM GPU 
    cudaMemcpy(h_B, d_B, num_bytes, cudaMemcpyDeviceToHost); 

    // Verify result is correct 
    CUTBoolean res = cutComparei((int *)B, (int *)h_B, N*N); 
    printf("Test %s\n",(1 == res)?"Passed":"Failed"); 
    printf("Elapsed Time for Sequential: \t%.2f ms\n", elapsed_time0); 
    printf("Elapsed Time for CUDA:\t%.2f ms\n", elapsed_time); 
    printf("CUDA Speedup:\t%.2fx\n",(elapsed_time0/elapsed_time)); 

    cudaFree(d_A); 
    cudaFree(d_B); 
    free(h_B); 

    cutilDeviceReset(); 
} 

Para la versión de memoria compartida, esto es lo que he probado hasta ahora:

#define N 1024 

__global__ void kernel(int *d_A, int *d_B, int width) { 
    //assuming width is 64 because that's the biggest number I can make it 
    //each MP has 48KB of shared mem, which is 12K ints, 32 threads/warp, so max 375 ints/thread? 
    __shared__ int A_sh[3][66]; 

    //get x and y index and turn it into linear index 

    for(i=0; i < width+2; i++) //have to load 2 extra values due to the -1 and +1 in algo 
      A_sh[index_y%3][i] = d_A[index+i-1]; //so A_sh[index_y%3][0] is actually d_A[index-1] 

    __syncthreads(); //and hope that previous and next row have been loaded by other threads in the block? 

    //ignore boundary conditions because it's pseudocode 
    for(i=0; i < width; i++) 
     d_B[index+i] = A_sh[index_y%3][i] + A_sh[index_y%3][i+2] + A_sh[index_y%3-1][i+1] + A_sh[index_y%3+1][i+1]; 

} 

main(){ 
    //same init as above until threads/grid init 

    dim3 threadsperblk(32,16); 
    dim3 numblks(32,64); 

    kernel<<<numblks,threadsperblk>>>(d_A,d_B,64); 

    //rest is the same 
} 

accidentes de este código mem compartido ("lanzamiento falló debido a la error no especificado ") ya que aún no he entendido todas las condiciones de contorno, pero no me preocupa tanto como encontrar la forma correcta de poner en marcha las cosas. Siento que mi código es demasiado complicado para ser el camino correcto (especialmente comparado con los ejemplos de SDK), pero tampoco veo otra forma de hacerlo, ya que mi matriz no encaja en las memorias compartidas como todos los ejemplos que puedo encontrar.

Y francamente, no estoy seguro de que sería mucho más rápido que en mi hardware (GTX 560 Ti - ejecuta la versión de la memoria global en 0.121ms), pero tengo que demostrar a mí mismo en primer lugar: P

Edición 2: para cualquier persona que se encuentre con esto en el futuro, el código en la respuesta es un buen punto de partida si desea hacer algo de memoria compartida.

Answer 1

La clave para sacar el máximo provecho de este tipo de operadores de esténciles en CUDA es la reutilización de datos. He descubierto que el mejor enfoque es generalmente hacer que cada bloque "camine" a través de una dimensión de la cuadrícula. Después de que el bloque haya cargado un mosaico inicial de datos en la memoria compartida, solo debe leerse una dimensión única (por lo tanto, una fila en un problema 2D de orden mayor de fila) de la memoria global para tener los datos necesarios en la memoria compartida para el segundo y posterior cálculos de filas El resto de los datos solo pueden ser reutilizados.Para visualizar cómo el buffer de memoria compartida se ve a través de los cuatro primeros pasos de este tipo de algoritmo:

"filas"

1. Tres (a, b, c) de la rejilla de entrada se cargan a la memoria compartida, y la plantilla calculado para fila (b) y se escribe en la memoria global

   AAAAAAAAAAAAAAAA bbbbbbbbbbbbbbbb cccccccccccccccc

2. Otra fila (d) se carga en la memoria intermedia de la memoria compartida, en sustitución de la fila (a), y los cálculos realizados para la fila (c) usando una galería de símbolos diferente, reflejando dónde los datos de la fila i s en la memoria compartida

   dddddddddddddddd bbbbbbbbbbbbbbbb cccccccccccccccc

3. Otra fila (e) se carga en la memoria intermedia de la memoria compartida, en sustitución de la fila (B), y los cálculos realizados para la fila (d), utilizando un diferente stencil de cualquiera de la etapa 1 o 2.

   dddddddddddddddd eeeeeeeeeeeeeeee cccccccccccccccc

4. Otra fila (f) se carga en la memoria intermedia de la memoria compartida, reemplazando la fila (c), y los cálculos hechos para la fila (e). Ahora los datos vuelven al mismo diseño que se utilizó en el paso 1, y se puede usar el mismo stencil utilizado en el paso 1.

   dddddddddddddddd eeeeeeeeeeeeeeee ffffffffffffffff

todo el ciclo se repite hasta que el bloque tiene atravesar longitud de la columna completa de la red de entrada. La razón para usar plantillas diferentes en lugar de cambiar los datos en el búfer de memoria compartida se reduce al rendimiento: la memoria compartida solo tiene un ancho de banda de 1000 Gb/s en Fermi, y el cambio de datos se convertirá en un cuello de botella en un código totalmente óptimo. Debe probar diferentes tamaños de búfer, ya que puede encontrar búferes más pequeños que permiten una mayor ocupación y un mejor rendimiento del kernel.

EDIT: Para dar un ejemplo concreto de cómo podría implementarse:

template<int width> 
__device__ void rowfetch(int *in, int *out, int col) 
{ 
    *out = *in; 
    if (col == 1) *(out-1) = *(in-1); 
    if (col == width) *out(+1) = *(in+1); 
} 

template<int width> 
__global__ operator(int *in, int *out, int nrows, unsigned int lda) 
{ 
    // shared buffer holds three rows x (width+2) cols(threads) 
    __shared__ volatile int buffer [3][2+width]; 

    int colid = threadIdx.x + blockIdx.x * blockDim.x; 
    int tid = threadIdx.x + 1; 

    int * rowpos = &in[colid], * outpos = &out[colid]; 

    // load the first three rows (compiler will unroll loop) 
    for(int i=0; i<3; i++, rowpos+=lda) { 
     rowfetch<width>(rowpos, &buffer[i][tid], tid); 
    } 

    __syncthreads(); // shared memory loaded and all threads ready 

    int brow = 0; // brow is the next buffer row to load data onto 
    for(int i=0; i<nrows; i++, rowpos+=lda, outpos+=lda) { 

     // Do stencil calculations - use the value of brow to determine which 
     // stencil to use 
     result =(); 
     // write result to outpos 
     *outpos = result; 

     // Fetch another row 
     __syncthreads(); // Wait until all threads are done calculating 
     rowfetch<width>(rowpos, &buffer[brow][tid], tid); 
     brow = (brow < 2) ? (brow+1) : 0; // Increment or roll brow over 
     __syncthreads(); // Wait until all threads have updated the buffer 
    } 
}