CUDA：在一个相对大的2D数组上使用共享内存

Question

CUDA：在一个相对大的2D数组上使用共享内存

cudagpu-shared-memory

5

我在课程作业中遇到了一个简单的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
}

这个共享内存代码会崩溃（“由于未指定的错误而启动失败”），因为我还没有捕获所有边界条件，但我不太担心这个问题，而是要找到正确的方法来让事情顺利进行。我觉得我的代码过于复杂，不应该是正确的路径（特别是与SDK示例相比），但我也看不到其他的方法，因为我的数组不像我能找到的所有示例那样适合放入共享内存中。

老实说，在我的硬件上（GTX 560 Ti - 全局内存版本运行时间为0.121ms），我不确定它会更快，但我需要先证明给自己看:P

编辑2：对于将来遇到此问题的任何人，答案中的代码是一个很好的起点，如果您想使用一些共享内存。

- a5ehren

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- talonmies · Accepted Answer

在CUDA中，获取最大性能的关键是数据的重复利用。通常最好的做法是让每个块“遍历”一维网格。在将初始数据块加载到共享内存后，只需要从全局内存中读取单个维度（例如，在按行优先顺序处理2D问题时为行）即可为第二行及以后的计算提供所需的共享内存中的数据。其余数据可以被重复利用。以下是通过这种算法的前四步可视化共享内存缓冲区的方式：

1. 将输入网格的三个“行”（a、b、c）加载到共享内存中，并计算出行(b)的模板并写入全局内存。

aaaaaaaaaaaaaaaa bbbbbbbbbbbbbbbb cccccccccccccccc

2. 另一行(d)被加载到共享内存缓冲区中，代替行(a)，并使用不同的模板为行(c)进行计算，反映了行数据在共享内存中的位置。

dddddddddddddddd bbbbbbbbbbbbbbbb cccccccccccccccc

3. 另一行(e)被加载到共享内存缓冲区中，代替行(b)，并使用步骤1或步骤2中的不同模板对行(d)进行计算。

dddddddddddddddd eeeeeeeeeeeeeeee cccccccccccccccc

4. 另一行(f)被加载到共享内存缓冲区中，代替行(c)，并使用不同的模板对行(e)进行计算。现在数据又回到了与步骤1中使用的相同布局，并且可以使用步骤1中使用的相同模板。

dddddddddddddddd eeeeeeeeeeeeeeee ffffffffffffffff

整个循环会重复进行，直到块遍历了输入网格的完整列长度。之所以使用不同的模板而不是将数据移动到共享内存缓冲区中，是因为性能原因-在Fermi上，共享内存只有约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
    }
}