Пример программы

Рассмотрим пример задачи, которая очень хорошо параллелится. Это может быть перемножение матриц или сложение векторов (об этом можно прочитать почти во всех статьях про CUDA). А мы посчитаем интеграл 4/(1+x*x) на отрезке [0,1] методом прямоугольников. Как нам подсказывают, он равен Пи. Вот так может выглядеть простейший код на С:

int main() { int numSteps = 10000; double left = 0.0; double right = 1.0; double step = (right-left)/numSteps; double sum = 0; for (double x = left + 0.5*step; x < right; x += step) sum += 4.0/(1.0 + x*x); printf("%0.8f.\n", sum/numSteps); return 0; }

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

Интегрирование по формуле прямоугольников: каждая полоска соответствует "шагу" фиксированной ширины. Высота каждой полоски равна значению подынтегральной функции. Если соединить вместе все полоски, можно приблизительно вычислить площадь под кривой, то есть значение интеграла. Создадим массив длиной numSteps и в каждый элемент этого массива запишем площадь соответствующего прямоугольника. Мы можем вычислить площадь одного прямоугольника независимо от других. Вот здесь нам и понадобится параллелизм. За вычисление площади каждого прямоугольника будет отвечать один поток. Далее нам остается только просуммировать элементы массива и получить значение интеграла.

double *a_h; // указатель на область памяти хоста double *a_d; // указатель на область памяти устройства a_h = (double *)malloc(sizeof(double)*numSteps); cudaMalloc((void **) &a_d, sizeof(double)*numSteps); int blockSize = 4; int blocks = numSteps / blockSize + (numSteps % blockSize == 0 ? 0:1);

cudaMalloc() - аналог сишной malloc() и она занимается тем, что выделяет область в глобальной памяти. Размер блока (blockSize) лучше всего делать равным 2^n - 4 , а структуры лучше выравнивать по границе в 16 байт (чтобы избежать так называемых конфликтов банков памяти).

Далее подготовим данные - занесем значения в массив.

double left = 0.0; double right = 1.0; double step = (right-left)/numSteps; int i = 0; for (double x = left + 0.5*step; x < right; x += step) { a_h[i] = x; i++; }

И передадим этот массив в память видеокарты.

cudaMemcpy(a_d, a_h, sizeof(double)*numSteps, cudaMemcpyHostToDevice);

Основные вычисления:

__global__ void calc(double *a, int n) { int idx = blockIdx.x * blockDim.x + threadIdx.x; double val = a[idx]; if (idx < n){ a[idx] = 4.0 /(1.0 + val*val); } }

Спецификатор __global__ показывает, что функция относится к ядру - её вызовет CPU, а выполнит GPU. Так же есть __device__ функция, которая выполнится на GPU и вызвать её можно только с GPU. Можно еще писать (а можно и не писать) спецификатор __host__ - функция вызывается CPU и выполняется на CPU, т.е. это - обычная функция. __global__ и __device__ функции не могут быть рекурсивными и должны содержать постоянное число аргументов. Т.к. функции __global__ и __device__ выполняются на GPU, то запустить их под обычным отладчиком и получить их адреса не получится. У NVIDIA есть специальные средства для этого, можно посмотреть на официальном сайте.
Каждый вызов __global__ функции должен соответсвовать спецификации вызова. Спецификация определяет размерность сетки и блоков, которые будут использоваться для выполнения этой функции на устройстве. Вызов должен соответсвовать форме:
     func<<< Dg, Db, Ns, S >>>(arguments)
Dg имеет тип dim3 и определяет размерность сетки, так Dg.x * Dg.y равно числу блоков. Тип dim3 трехмерный, но координата Dg.z обычно не используется.
Db тоже имеет тип dim3 и означает размерность и размер каждого блока. Значение Db.x * Db.y * Db.z равно числу потоков в блоке.
Ns имеет тип size_t и определяет число байтов в shared памяти, которая динамически размещается для каждого блока в дополнение к статической памяти. Ns необязательный параметр и по умолчанию равен 0.
Параметр S типа cudaStream_t , определяющий дочерние потоки. S также необязателен с параметром по умолчанию, равным нулю.
Встроенные переменные:
      blockIdx - номер блока внутри сетки
      threadIdx - номер потока внутри блока
      blockDim - число потоков в блоке
      blockIdx и blockDim- трехмерны и содержат поля x,y,z, а сама сетка двумерна.
Т.к. массивы у нас одномерные, используется только координата x. После того, как провели вычисления, нужно передать данные обратно на хост :
cudaMemcpy(a_h, a_d, size, cudaMemcpyDeviceToHost).
cudaMemcpyDeviceToHost - копирование с устройства на хост, cudaMemcpyHostToDevice - соответственно обратно. Еще несколько действий и выводим результат : 3,141592.
Полный текст программы можно посмотреть тут, откомпилированную версию программы здесь. Так же доступны для скачивания программа для устройств CUDA и её аналог в режиме эмуляции.