Блог

У меня была простая проблема с 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) amp;amp; (index_y > 0) amp;amp; (index_x < N-1) amp;amp; (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("%dn",A[i][j]);
        }
    }

    cudaEvent_t start_event0, stop_event0;
    float elapsed_time0;
    CUDA_SAFE_CALL( cudaEventCreate(amp;start_event0) );
    CUDA_SAFE_CALL( cudaEventCreate(amp;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(amp;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**)amp;d_A, num_bytes);
    cudaMalloc((void**)amp;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(amp;start_event) );
    CUDA_SAFE_CALL( cudaEventCreate(amp;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(amp;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 %sn",(1 == res)?"Passed":"Failed");
    printf("Elapsed Time for Sequential: t%.2f msn", elapsed_time0);
    printf("Elapsed Time for CUDA:t%.2f msn", elapsed_time);
    printf("CUDA Speedup:t%.2fxn",(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), но я также не вижу другого способа сделать это, поскольку мой массив не помещается в разделяемую память, как и все примеры, которые я могу найти.

И, честно говоря, я не уверен, что это было бы намного быстрее на моем оборудовании (GTX 560 Ti — запускает версию глобальной памяти за 0,121 мс), но сначала мне нужно доказать это самому себе: P

Редактировать 2: Для любого, кто столкнется с этим в будущем, код в ответе является хорошей отправной точкой, если вы хотите использовать некоторую общую память.

Ключом к получению максимальной отдачи от такого рода трафаретных операторов в CUDA является повторное использование данных. Я обнаружил, что наилучший подход обычно заключается в том, чтобы каждый блок «проходил» через измерение сетки. После того, как блок загрузил начальный фрагмент данных в разделяемую память, из глобальной памяти необходимо прочитать только одно измерение (строка в строке — проблема 2D основного порядка), чтобы иметь необходимые данные в разделяемой памяти для вычислений второй и последующих строк. Остальные данные можно просто использовать повторно. Чтобы визуализировать, как выглядит буфер разделяемой памяти, выполните первые четыре шага такого алгоритма:

1. Три «строки» (a, b, c) входной сетки загружаются в разделяемую память, а шаблон вычисляется для строки (b) и записывается в глобальную память

   ааааааааааааааааа ббббббббббббббббк cccccccccccccccccccccccccccc

2. В буфер общей памяти загружается другая строка (d), заменяющая строку (a), а вычисления, выполненные для строки (c), выполняются с использованием другого шаблона, отражающего, где данные строки находятся в общей памяти

   дддддддддддддддд ббббббббббббббббк cccccccccccccccccccccccccccccc

3. В буфер разделяемой памяти загружается другая строка (e), заменяющая строку (b) и вычисления, выполненные для строки (d), с использованием шаблона, отличного от шага 1 или 2.

   dddddddddddddddd eeeeeeeeeeeeeeee cccccccccccccccc

4. В буфер общей памяти загружается другая строка (f), заменяющая строку (c) и вычисления, выполненные для строки (e). Теперь данные возвращаются к тому же макету, который использовался на шаге 1, и можно использовать тот же шаблон, который использовался на шаге 1.

   dddddddddddddddd eeeeeeeeeeeeeeee ffffffffffffffff

Весь цикл повторяется до тех пор, пока блок не пройдет всю длину столбца входной сетки. Причина использования разных шаблонов вместо переноса данных в буфере разделяемой памяти связана с производительностью — пропускная способность разделяемой памяти в Fermi составляет всего около 1000 Гбит / с, и перенос данных станет узким местом в полностью оптимальном коде. Вам следует попробовать буферы разных размеров, потому что вы можете обнаружить, что буферы меньшего размера обеспечивают более высокую занятость и улучшенную пропускную способность ядра.

РЕДАКТИРОВАТЬ: чтобы привести конкретный пример того, как это может быть реализовано:

 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 = amp;in[colid], * outpos = amp;out[colid];

    // load the first three rows (compiler will unroll loop)
    for(int i=0; i<3; i  , rowpos =lda) {
        rowfetch<width>(rowpos, amp;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 = resu<

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