CUDA: mémoire partagée sur un tableau 2D de grande taille

Question

J'ai eu un problème CUDA simple pour une affectation de classe, mais le professeur a ajouté une tâche facultative pour implémenter le même algorithme en utilisant la mémoire partagée à la place. Je ne pouvais pas le finir avant la date limite (comme dans, la date d'entrée était il y a une semaine) mais je suis toujours curieux alors maintenant je vais demander à Internet;). La tâche de base consistait à implémenter une version bastardisée d'une sur-relaxation successive rouge-noire à la fois séquentiellement et dans CUDA, assurez-vous d'avoir le même résultat dans les deux et de comparer ensuite l'accélération. Comme je l'ai dit, le faire avec de la mémoire partagée était un ajout facultatif de + 10%. Je vais poster ma version de travail et pseudocode ce que j'ai essayé de faire puisque je n'ai pas le code dans mes mains pour le moment, mais je peux le mettre à jour plus tard avec le code réel si quelqu'un a besoin il.

Avant que quelqu'un ne le dise: Oui, je sais que CUtil est boiteux, mais cela facilite la comparaison et les temporisations.

version mémoire globale de travail:

#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(); 
} 

Pour la version de mémoire partagée, ce que j'ai essayé jusqu'à présent:

#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 
} 

Ce plantages de code mem partagé ("lancement a échoué en raison de erreur non spécifiée ") puisque je n'ai pas encore pris toutes les conditions aux limites, mais je ne m'inquiète pas autant que de trouver la bonne façon de faire avancer les choses. Je pense que mon code est trop compliqué pour être le bon chemin (surtout par rapport aux exemples du SDK), mais je ne peux pas non plus voir une autre façon de le faire puisque mon tableau ne rentre pas dans mem comme tous les exemples peuvent trouver.

Et franchement, je ne suis pas sûr que ce serait beaucoup plus rapide sur mon matériel (GTX 560 Ti - exécute la version de la mémoire globale en 0.121ms), mais je dois me prouver d'abord: P

Édition 2: Pour ceux qui courent à travers cela à l'avenir, le code dans la réponse est un bon point de départ si vous voulez faire de la mémoire partagée.

Answer 1

La réutilisation des données est la clé pour tirer le meilleur parti de ces types d'opérateurs de stencil dans CUDA. J'ai trouvé que la meilleure approche consiste habituellement à faire «marcher» chaque bloc à travers une dimension de la grille. Après que le bloc a chargé une mosaïque de données initiale dans la mémoire partagée, une seule dimension (donc rangée dans une rangée - problème majeur de l'ordre 2D) doit être lue dans la mémoire globale pour avoir les données nécessaires dans la mémoire partagée calculs de rangées. Le reste des données peut simplement être réutilisé.Pour visualiser la façon dont la mémoire tampon de la mémoire partagée regarde à travers les quatre premières étapes de ce genre d'algorithme:

1. Trois « rangées » (a, b, c) de la grille d'entrée sont chargés dans la mémoire partagée, et le pochoir calculé pour la ligne (b) et écrites dans la mémoire globale

   aaaaaaaaaaaaaaaa bbbbbbbbbbbbbbbb cccccccccccccccc

2. autre rangée (d) est chargé dans la mémoire tampon de la mémoire partagée, en remplacement de la ligne (a), et les calculs effectués pour la rangée (c) en utilisant un pochoir différent, reflétant où les données de ligne i s dans la mémoire partagée

   DDDDDDDDDDDDDDDD bbbbbbbbbbbbbbbb cccccccccccccccc

3. autre rangée (e) est chargé dans la mémoire tampon de la mémoire partagée, en remplacement de la ligne (b), et les calculs effectués pour la ligne (d), en utilisant un autre stencil soit de l'étape 1 ou 2.

   DDDDDDDDDDDDDDDD eeeeeeeeeeeeeeee cccccccccccccccc

4. autre rangée (f) est chargé dans la mémoire tampon de la mémoire partagée, remplacer la rangée (c) et les calculs effectués pour la rangée (e). Maintenant, les données sont de retour à la même mise en page que celle utilisée à l'étape 1, et le même stencil utilisé à l'étape 1 peut être utilisé.

   DDDDDDDDDDDDDDDD eeeeeeeeeeeeeeee ffffffffffffffff

L'ensemble du cycle se répète jusqu'à ce que le bloc a traverser toute la longueur de la colonne de la grille d'entrée. La raison d'utiliser différents stencils plutôt que de déplacer les données dans le tampon de mémoire partagée est la performance - la mémoire partagée n'a qu'une bande passante d'environ 1000 Gb/s sur Fermi, et le décalage des données devient un goulot d'étranglement dans le code optimal. Vous devriez essayer différentes tailles de tampons, car vous pourriez trouver des tampons plus petits permettant une occupation plus élevée et un débit de noyau amélioré.

EDIT: Pour donner un exemple concret de la façon dont cela pourrait être mis en œuvre:

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 
    } 
}