CUDA: общая память через двумерный массив большого размера

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
}

Эта общая памятьпроисходит сбой кода («запуск не удался из-за неуказанной ошибки»), поскольку я еще не уловил все граничные условия, но я не столько беспокоюсь об этом, сколько ищу правильный способ развития событий.Я чувствую, что мой код слишком сложен, чтобы быть правильным путем (особенно по сравнению с примерами SDK), но я также не вижу другого способа сделать это, так как мой массив не вписывается в общую память, как все примеры, которые яможно найти.

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

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

Answer 1

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

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

   aaaaaaaaaaaaaaaa bbbbbbbbbbbbbbbb cccccccccccccccc

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

   dddddddddddddddd bbbbbbbbbbbbbbbb cccccccccccccccc

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

   dddddddddddddddd eeeeeeeeeeeeeeee cccccccccccccccc

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

   dddddddddddddddd eeeeeeeeeeeeeeee ffffffffffffffff

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

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

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