CUDA：通過大型2D陣列共享內存

Question

我對類賦值有一個簡單的CUDA問題，但教授添加了一個可選任務來實現使用共享內存的相同算法。我無法在截止日期之前完成它（如轉機日期是一週前），但我仍然好奇，所以現在我要問網絡了）。

基本的任務是實現連續過度放鬆和CUDA過程中的混合版本的連續過度放鬆，確保你在兩者中得到相同的結果，然後比較加速比。就像我說的那樣，使用共享內存做一個可選的+ 10％附加組件。

我打算髮布我的工作版本和僞代碼，因爲我目前沒有手中的代碼，但如果有人需要，我可以稍後用實際代碼更新此代碼它。

之前有人說：是的，我知道使用CUtil是跛腳，但它使比較和計時器更容易。

工作全局存儲器版本：

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

共享內存的版本，這是我到目前爲止已經試過：

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

此共享MEM代碼崩潰（「因發射失敗未指定的錯誤「），因爲我還沒有捕捉到所有的邊界條件，但我並不擔心這一點，而是尋找正確的方式來推動事情發展。我覺得我的代碼太複雜了，不能成爲正確的路徑（特別是與SDK例子相比），但是我也看不到另一種方法來實現它，因爲我的數組不像所有例子那樣適合共享mem可以找到。

坦率地說，我不知道這將是我的硬件要快得多（GTX 560 TI - 運行在0.121ms全局內存版本），但我需要首先證明自己：P

編輯2：對於任何在將來運行的人來說，如果你想做一些共享內存，答案中的代碼是一個很好的起點。

Answer 1

在CUDA中獲得這些模板操作符的最大好處的關鍵是重用數據。我發現最好的方法通常是讓每個塊在網格的一個維度上「走」。塊將數據的初始數據塊加載到共享內存之後，需要從全局內存中僅讀取一個維度（行中的行爲主要次序2D問題），以便在第二個及隨後的共享內存中擁有必要的數據行計算。剩下的數據可以被重用。爲了顯現共享存儲器緩衝區的外觀通過所述前四個步驟這種算法的：

1. 三「行」（A，B，C）的輸入網格的被加載到共享存儲器中，並且所述模板計算用於行（b）和寫入到全局存儲器

   aaaaaaaaaaaaaaaa bbbbbbbbbbbbbbbb cccccccccccccccc

2. 另一行（d）被加載到共享存儲器緩衝區，替換行（a）中，並且計算用於行由（c）使用不同的模板，反映行數據i S IN共享存儲器

   dddddddddddddddd bbbbbbbbbbbbbbbb cccccccccccccccc

3. 另一行（e）中被加載到共享存儲器緩衝區，替換行（b）和用於行（d）進行的計算，使用不同的從步驟1或2

   模版dddddddddddddddd eeeeeeeeeeeeeeee cccccccccccccccc

4. 另一行（F）被加載到共享存儲器緩衝區，替換行（c）和對行（e）進行的計算。現在數據返回到步驟1中使用的相同佈局，並且可以使用步驟1中使用的相同模板。

   dddddddddddddddd eeeeeeeeeeeeeeee FFFFFFFFFFFFFFFF

整個週期反覆進行，直到塊具有穿過所述輸入格的整列長度。之所以使用不同的模具，而不是在所述共享存儲器緩衝區中的數據移位下降到性能 - 共享存儲器上費米只有約1000 Gb/s的帶寬，和數據的移動將成爲充分最佳代碼的瓶頸。您應該嘗試不同的緩衝區大小，因爲您可能會發現較小的緩衝區允許更高的佔用率和更高的內核吞吐量。

編輯：爲了舉例說明如何可能實現一個具體的例子：

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