Простейшие примеры использования CUDA

Что такое CUDA? Смотрите в википедии.
Рассмотрим коротко простейшие примеры использования этой интересной технологии. Этот пост является конспектом презентации Introduction to CUDA C автора Jason Sanders на конференции NVIDIA GTC 2010 - GPU Technology Conference. (Pre-Conference Tutorials, GTC 2010 - ID 2131).
Итак. Самая базовая терминология:
Host - CPU и его память (host memory);
Device - GPU и его память (device memory);
Начнем с “Hello, World!”. А куда ж без него?
Код делится на тот, что выполняется на CPU и тот что выполняется на GPU. Следующий самый просто пример, это программа, которая выводит приветствие миру и вызывает функцию из GPU, которая ничего не делает:

__global__ void kernel( void ) {
}

int main( void ) {
        kernel<<< 1, 1 >>>();
        printf( "Hello, World!\n" );
        return 0;
}

Ключевое слово __global__ в объявлении и определении функции обозначает, что функция выполняется на графической карте. И что эта функция вызывается из CPU. Неизменным аттрибутом такой функции являются скобки <<< >>> при вызове этой функции. Аргументы внутри этих скобок обсудим ниже. Компилятор nvcc разделяет исходный код на части, которые выполняются отдельно на хосте и девайсе. Код, сгенерированный для таких функций как kernel из данного примера, выполняется на девайсе. Кстати, как я понял, device - это не обязательно GPU.
Рассмотрим более сложный пример. CPU передает некие величины в GPU для сложения их там.

__global__ void add( int *a, int *b, int *c ) {
    *c = *a + *b;
}

int main( void ) {
    int a, b, c;                   // host копии a, b, c
    int *dev_a, *dev_b, *dev_c;    // device копии of a, b, c
    int size = sizeof( int );      
    //выделяем память для device копий для a, b, c
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, size );
    a = 2;
    b = 7;
    // копируем ввод на device
    cudaMemcpy( dev_a, &a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, &b, size, cudaMemcpyHostToDevice );
    // запускаем add() kernel на GPU, передавая параметры
    add<<< 1, 1 >>>( dev_a, dev_b, dev_c );
    // copy device result back to host copy of c
    cudaMemcpy( &c, dev_c, size, cudaMemcpyDeviceToHost );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

Базовый CUDA API для работы с памятью GPU (и девайса вообще):
cudaMalloc();
cudaFree();
cudaMemcpy();
Их Си аналоги: malloc(), free(), memcpy().

__global__ void add( int *a, int *b, int *c ) {
    c[blockIdx.x] = a[blockIdx.x] + b[blockIdx.x];
}

#define N  512
int main( void ) {
    int *a, *b, *c;                   // host копии a, b, c
    int *dev_a, *dev_b, *dev_c;       // device копии a, b, c
    int size = N * sizeof( int );     
    //выделяем память для device копий a, b, c
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, size );
    a = (int*)malloc( size );
    b = (int*)malloc( size );
    c = (int*)malloc( size );
    random_ints( a, N );
    random_ints( b, N );
    // копируем ввод на device
    cudaMemcpy( dev_a, a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, size, cudaMemcpyHostToDevice );
    // launch add() kernel with N parallel blocks
    add<<< N, 1 >>>( dev_a, dev_b, dev_c );
    // копируем результат работы device обратно на host - копию c
    cudaMemcpy( c, dev_c, size, cudaMemcpyDeviceToHost );
    free( a ); free( b ); free( c );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

Главное отличие этого кода от предыдущих это аргумент N в вызове функции add<<< N, 1 >>>. Он обозначает число так называемых блоков, на которые разбивается код на девайсе при распараллеливании. Обращаться к отдельному блоку можно через индекс массива блоков blockIdx.x (это встроенная переменная). Блоки, в свою очередь могут разбиваться на нити (threads)
Перепишем наш код для использования нитей, а не блоков (будем использовать не N блоков по одной нити, а один блок с N нитями):

__global__ void add( int *a, int *b, int *c ) {
     c[ threadIdx.x ] = a[ threadIdx.x ] + b[ threadIdx.x ];
     //    blockIdx.x         blockIdx.x         blockIdx.x
}

#define N  512
int main( void ) {
    int *a, *b, *c;                      //host копии a, b, c
    int *dev_a, *dev_b, *dev_c;          //device копии of a, b, c
    int size = N * sizeof( int );        
    //выделяем память для копий a, b, c
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, size );
    a = (int*)malloc( size );
    b = (int*)malloc( size );
    c = (int*)malloc( size );
    random_ints( a, N );
    random_ints( b, N );
    // копируем ввод device
    cudaMemcpy( dev_a, a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, size, cudaMemcpyHostToDevice );
    // запускаем на выполнение add() kernel с N нитями в блоке
    add<<< 1, N >>>( dev_a, dev_b, dev_c );
    //     N, 1
    // копируем результат работы device обратно на host (копия c)
    cudaMemcpy( c, dev_c, size, cudaMemcpyDeviceToHost );
    free( a ); free( b );  free( c );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

Теперь рассмотрим код, который использует и несколько блоков и несколько нитей в каждом блоке:

__global__ void add( int *a, int *b, int *c ) {
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    c[index] = a[index] + b[index];
}

#define N  (2048*2048)
#define THREADS_PER_BLOCK 512
int main( void ) {
    int *a, *b, *c;                      // host копии a, b, c
    int *dev_a, *dev_b, *dev_c;          // device копии of a, b, c
    int size = N * sizeof( int );        
    //выделяем память на device для of a, b, c
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, size );
    a = (int*)malloc( size );
    b = (int*)malloc( size );
    c = (int*)malloc( size );
    random_ints( a, N );
    random_ints( b, N );
    //копируем ввод на device
    cudaMemcpy( dev_a, a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, size, cudaMemcpyHostToDevice );
    //запускаем на выполнение add() kernel с блоками и нитями
    add<<< N/THREADS_PER_BLOCK, THREADS_PER_BLOCK >>>( dev_a, dev_b, dev_c );
    // копируем результат работы device на host ( копия c )
    cudaMemcpy( c, dev_c, size, cudaMemcpyDeviceToHost );
    free( a ); free( b ); free( c );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

Нити на первый взгляд кажутся избыточными. Параллельные нити, в отличие параллельных блоков, имеют механизмы для, так называемых, коммуникации и синхронизации. Давайте рассмотрим этот вопрос.
Возьмем в качестве примера на этот раз не сложение, а скалярное произведение.
Нити в одном блоке могут иметь общую память shared memory. Эта память декларируется в коде ключевым словом __shared__. Эта память не видна в другом блоке.
Также нити можно "синхронизировать". То есть каждая нить в некой "точке синхронизации" будет ждать пока выполнение других нитей тоже не достигнет этой точки. Синхронизация нитей истинна только для нитей одного блока. Синхронизация вводиться функцией __syncthreads().
Вот код, с синхронизацией и использованием shared memory:

//__global__ void dot( int *a, int *b, int *c )   {
    // Каждая нить вычисляет одно слагаемое скалярного произведения
//    int temp = a[threadIdx.x] * b[threadIdx.x];
    //Не может вычислить сумму
    //Каждая копия ‘temp’ приватна для своей нити
//}

__global__ void dot( int *a, int *b, int *c ) {
      __shared__ int temp[N];
     temp[threadIdx.x] = a[threadIdx.x] * b[threadIdx.x];
      __syncthreads();
      if( 0 == threadIdx.x ) {
          int sum = 0;
          for( int i = 0; i < N; i++ )
              sum += temp[i];
          *c = sum;
      }
}

#define N  512
int main( void ) {
    int *a, *b, *c;                      
    int *dev_a, *dev_b, *dev_c;          
    int size = N * sizeof( int );        
    
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, sizeof( int ) );
    a = (int *)malloc( size );
    b = (int *)malloc( size );
    c = (int *)malloc( sizeof( int ) );
    random_ints( a, N );
    random_ints( b, N );
    // копируем ввод на device
    cudaMemcpy( dev_a, a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, size, cudaMemcpyHostToDevice );
    //запускаем на выполнение dot() kernel с 1 блоком и N нитями
    dot<<< 1, N >>>( dev_a, dev_b, dev_c );
    //копируем результат работы device на host копией c
    cudaMemcpy( c, dev_c, sizeof( int ) , cudaMemcpyDeviceToHost );
    free( a ); free( b ); free( c );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

threadIdx.x используется для доступе к отдельной нити. threadIdx - это встроенная переменная. Для индексации вводимых/выводимых элементов используем выражение (threadIdx.x + blockIdx.x * blockDim.x).blockDim.x - это встроенная переменная, обозначающая число нитей в блоке.
В последнем примере использовался только один блок. А если нам надо использовать много блоков? Скалярное произведение же у нас тут только для примера.

__global__ void dot( int *a, int *b, int *c ) {
    __shared__ int temp[THREADS_PER_BLOCK];
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    temp[threadIdx.x] = a[index] * b[index];
    __syncthreads();
    if( 0 == threadIdx.x ) {
        int sum = 0;
        for( int i = 0; i < THREADS_PER_BLOCK; i++ )
            sum += temp[i];
        *c += sum; 
    }
}

В таком варианте функции dot у нас возникает race condition. Что это? Не тема этого поста. Смотрим википедию. С этой проблемой боремся с помощью атомарных функций. В данном случае - это atomicAdd.

__global__ void dot( int *a, int *b, int *c ) {
     __shared__ int temp[THREADS_PER_BLOCK];
     int index = threadIdx.x + blockIdx.x * blockDim.x;
     temp[threadIdx.x] = a[index] * b[index];
     __syncthreads();
     if( 0 == threadIdx.x ) {
         int sum = 0;
         for( int i = 0; i < THREADS_PER_BLOCK; i++ )
             sum += temp[i];
         atomicAdd( c , sum );
     }
}

Функция main будет такая:

#define N (2048*2048)
#define THREADS_PER_BLOCK 512
int main( void ) {
    int *a, *b, *c;                      // host копии a, b, c
    int *dev_a, *dev_b, *dev_c;          // device копии a, b, c
    int size = N * sizeof( int );        
    //выделяем место в памяти для device копий a, b, c
    cudaMalloc( (void**)&dev_a, size );
    cudaMalloc( (void**)&dev_b, size );
    cudaMalloc( (void**)&dev_c, sizeof( int ) );
    a = (int *)malloc( size );
    b = (int *)malloc( size );
    c = (int *)malloc( sizeof( int ) );
    random_ints( a, N );
    random_ints( b, N );
    //копируем ввод на device
    cudaMemcpy( dev_a, a, size, cudaMemcpyHostToDevice );
    cudaMemcpy( dev_b, b, size, cudaMemcpyHostToDevice );
    //запускаем на выполнение dot() kernel
    dot<<< N/THREADS_PER_BLOCK, THREADS_PER_BLOCK >>>( dev_a, dev_b, dev_c );
    //копируем результат работы device обратно на host (копия c)
    cudaMemcpy( c, dev_c, sizeof( int ) , cudaMemcpyDeviceToHost );
    free( a ); free( b ); free( c );
    cudaFree( dev_a );
    cudaFree( dev_b );
    cudaFree( dev_c );
    return 0;
}

CUDA предоставляет атомарные функции, которые гарантируют "безопасность" read-modify-write операции.
Пример атомарных функций в CUDA:
atomicInc()
atomicAdd()
atomicDec()
atomicSub()
atomicExch()
atomicMin()
atomicCAS()
atomicMax()
С помощью этих примеров продемонстрированы приемы, которые можно использовать для применения в параллельных вычислениях на видеокартах NVIDIA.