Основной принцип GPU-расчетов - параллельность вычислений. Данные, хранящиеся в “глобальной памяти” (global & constant memory) устройства, обрабатываются модулями (каждый модуль называется “ядром”), каждый из которых работает параллельно с другими. Модуль имеет и свою собственную память для промежуточных данных (private memory). Так это выглядит в виде блок-схемы:
Таким образом, если задача может быть разбита на небольшие блоки, параллельно обрабатывающие небольшой фрагмент блока данных, такая задача может эффективно быть решена на GPU.
Рассмотрим пример: необходимо проверить, какие числа в массиве являются простыми. Массив может быть большим, например миллион элементов. Такая задача идеальна для распараллеливания: каждое число может быть проверено независимо от предыдущего.
Для решения такой задачи с помощью OpenCL необходимо выполнить ряд шагов.
Написать код микроядра (kernel):
Этот код будет запускаться непосредственно на графических процессорах видеокарты. Код пишется на языке C. В данном примере мы для упрощения храним код прямо в виде строки в программе.
const char *KernelSource = "\n" \
"__kernel void primes( \n" \
" __global unsigned int* input, \n" \
" __global unsigned int* output) \n" \
"{ \n" \
" unsigned int i = get_global_id(0); \n" \
" //printf(\"Task-%d\\n\", i); \n" \
" output[i] = 0; \n" \
" unsigned int val = input[i]; \n" \
" for(unsigned int p=2; p<=val/2; p++) { \n" \
" if (val % p == 0) \n" \
" return; \n" \
" } \n" \
" output[i] = 1; \n" \
"} \n" \
"\n";
Суть кода проста. Массив input хранит числа, которые нужно проверить, функция get_global_id возвращает индекс задачи, которую выполняет данное ядро. Мы берем число с нужным индексом, проверяем его на простоту, и записываем 0 или 1 в зависимости от результата, в массив output.
Инициализировать подготовку вычислений:
int gpu = 1;
clGetDeviceIDs(NULL, gpu ? CL_DEVICE_TYPE_GPU : CL_DEVICE_TYPE_CPU, 1, &device_id, NULL);
cl_context context = clCreateContext(0, 1, &device_id, NULL, NULL, &err);
cl_command_queue commands = clCreateCommandQueue(context, device_id, 0, &err);
На этом этапе можно выбрать где будут производиться вычисления, на основном процессоре или на GPU. Для отладки удобнее основной процессор, окончательные расчеты быстрее на GPU.
Подготовить данные
#define DATA_SIZE 1024
cl_uint *data = (cl_uint*)malloc(sizeof(cl_uint) * DATA_SIZE);
cl_uint *results = (cl_uint*)malloc(sizeof(cl_uint) * DATA_SIZE);
Загрузить данные и программу из основной памяти в GPU
cl_program program = clCreateProgramWithSource(context, 1, (const char **) & KernelSource, NULL, &err);
clBuildProgram(program, 0, NULL, NULL, NULL, NULL);
cl_kernel kernel = clCreateKernel(program, "primes", &err);
cl_mem output = clCreateBuffer(context, CL_MEM_WRITE_ONLY, sizeof(cl_uint) * count, NULL, NULL);
clEnqueueWriteBuffer(commands, input, CL_TRUE, 0, sizeof(cl_uint) * count, data, 0, NULL, NULL);
clSetKernelArg(kernel, 0, sizeof(cl_mem), &output);
clGetKernelWorkGroupInfo(kernel, device_id, CL_KERNEL_WORK_GROUP_SIZE, sizeof(local), &local, NULL);
Запустить вычисления на GPU и дождаться их завершения
global = DATA_SIZE;
clEnqueueNDRangeKernel(commands, kernel, 1, NULL, &global, &local, 0, NULL, NULL);
clFinish(commands);
Загрузить результаты обратно из GPU в основную память
clEnqueueReadBuffer( commands, output, CL_TRUE, 0, sizeof(cl_uint) * count, results, 0, NULL, NULL );
Освободить данные
free(data);
free(results);
clReleaseMemObject(input);
clReleaseMemObject(output);
clReleaseProgram(program);
clReleaseKernel(kernel);
clReleaseCommandQueue(commands);
clReleaseContext(context);
Как можно видеть, процесс довольно-таки громоздкий, но оно того стоит. Для примера, проверка простоты 250000 чисел заняла на процессоре Core i5 около 6 секунд. И всего лишь 0.5 секунд заняло выполнение вышеприведенного кода на встроенной видеокарте. Для дешевого нетбука с процессором Intel Atom этот же код выполнялся 34 секунды на основном процессоре, и 6 секунд на GPU. Т.е. разница весьма прилична.
Разумеется, еще раз стоит повторить, что “игра стоит свеч” лишь в том случае, если задача хорошо распараллеливается на небольшие блоки, в таком случае выигрыш будет заметен.
Владельцы видеокарт NVIDIA (особенно игровых и достаточно мощных) могут также посмотреть в сторону библиотеки NVIDIA CUDA, расчеты с ее помощью должны быть еще быстрее.
Приложение 2 - Быстродействие языка Python
Язык Python очень удобен своей краткостью и лаконичностью, возможностью использования большого количества сторонних библиотек. Однако, один из его минусов, который может быть ключевым для математических расчетов - это быстродействие. Python это интерпретатор, он не создает exe-файл, что разумеется, сказывается на скорости выполнения программы.
Рассмотрим простой пример: рассчитаем сумму квадратов чисел от 1 до 1000000. Также выведем время выполнения программы.
Программа на языке Python выглядит так:
import time
start_time = time.time()
s = 0
for x in xrange(1,1000001):
s += x*x
print("Sum={}, T={}s".format(s, time.time() - start_time))
Результаты работы:
Sum = 333333833333500000, T = 0.47s
Учитывая, что чисел всего миллион, не так уж и быстро. Попробуем ускорить программу, для этого по возможности используем функции встроенных библиотек. Они зачастую написаны на С, и работают быстрее.
import time
start_time = time.time()
l = xrange(1000001)
s = sum(x*x for x in l)
print("Sum = {}, T = {}s".format(s, time.time() - start_time))
Результаты работы:
Sum = 333333833333500000, T = 0.32s
Быстрее, но лишь чуть-чуть.
И наконец, призываем “тяжелую артиллерию”: перепишем программу на языке Cи.
Код выглядит так:
#include <stdio.h>
#include <time.h>
int main()
{
clock_t start = clock();
unsigned long long int sum = 0, i;
for(i=1; i<1000001; i++) {
sum += i*i;
}
clock_t end = clock();
printf("Sum = %llu, T = %fs", sum, (float)(end - start)/CLOCKS_PER_SEC);
return 0;
}
Как можно видеть, он ненамного сложнее python-версии. Перед запуском программы, ее надо скомпилировать, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest". Результат очевиден: T = 0.007 секунд. И еще чуть-чуть: добавляем флаг оптимизации по скорости, выполнив команду C:\GCC\bin\gcc.exe "Appendix-2 - speedTest.c" -o"Appendix-2 - speedTest" -O3. Результат: 0.0035 секунд, разница в быстродействии более 100 раз!
Увы, в более сложных задачах такого прироста реально не бывает (в последнем примере очень короткий код, который видимо полностью помещается в кеш-памяти процессора), но на некоторое улучшение быстродействия можно рассчитывать. Хотя переписывание программы - это крайний случай, сначала целесообразно поискать стандартные библиотеки, которые возможно уже решают данную задачу. К примеру, следующий код на языке Python вычисляет сумму элементов массива за 0.1с:
a = xrange(1000001)
s = 0
for x in a:
s += x
print(s)
Можно использовать встроенную функцию sum:
a = range(1000001)
s = sum(a)
print(s)
Данный код выполняется за 0.02 секунды, т.е. в 5 раз быстрее первого варианта.
Вопрос, что выбрать, на самом деле, не так однозначен - язык Си быстрее, но и сложнее в использовании. Но если заранее известно, что задача состоит в обработке большого набора чисел (например поиск простых чисел или магических квадратов), то может быть более целесообразным сразу писать программу на Си или С++, в принципе это не намного сложнее, а работать программа будет быстрее. С другой стороны, при наличии в программе сложной логики или структур данных, написание ее на Python будет гораздо проще. Поэтому, на практике целесообразно прототип программы делать на языке Python, и только в том случае, если проблема быстродействия становится критичной и не решаемой доступными средствами, переходить на более низкоуровневые средства разработки.
Кстати, для тех, кто захочет максимально широко использовать Python в математических расчетах, крайне рекомендуется ознакомиться с библиотекой numpy, она является стандартом де-факто для всех серьезных расчетов.
Приложение 3 - Пример оптимизации программы
Рассмотрим практический пример сравнения скорости вычислений на Python. Допустим, есть массив целых чисел, для каждого числа надо найти количество его делителей.
Функцию нахождения делителей числа будем использовать “как есть”, без какой-либо оптимизации: перебираем все числа, если число делится без остатка, увеличиваем счетчик.
def get_dividers(n):
divs = 0
for i in xrange(1, n+1):
if n % i == 0:
divs += 1
return divs
Рассмотрим способы решения задачи. Исходные данные: массив целых чисел:
values = [41212317, 4672313, 4342311, 46512319, 51212317, 5672313, 5342311, 56512319]
Способ-1
Решение “в лоб” (заодно показан способ измерения времени выполнения).
import timeit
start_time = timeit.default_timer()
res = []
for v in values:
res.append(get_dividers(v))
print(res)
print("T =", timeit.default_timer() - start_time)
Время выполнения кода на компьютере с процессором Core i7: 10c.
Способ-2
Используем функцию map, позволяющую применить функцию сразу к массиву.
res = map(get_dividers, values)
print(list(res))
Время выполнения: те же 10с, пока мы выиграли разве что в краткости записи.
Способ-3
Используем многопоточность, класс multiprocessing.Pool, позволяющий разбить вычисления над массивом на несколько потоков.
p = multiprocessing.Pool(processes=4)
res = p.map(get_dividers, values)
print(list(res))
Время выполнения: 3.7с, что уже лучше. Интересно, что несмотря на разбивку на 4 процесса, реальное время выполнения уменьшилось лишь вдвое. Многопоточные вычисления - достаточно дорогостоящая в плане накладных расходов операция, множество ресурсов процессора тратится на синхронизацию и передачу данных между процессами.
Способ-4
Используем библиотеку параллельных вычислений ipyparallel. Она имеет большие возможности, даже позволяет производить вычисления параллельно на нескольких компьютерах.
Перед тестом запустим сервер вычислений на том же компьютере, командой
ipcluster start -n 4.
Код запуска вычислений не намного больше:
import ipyparallel as ipp
rc = ipp.Client()
dv = rc[:]
res = dv.map_sync(get_dividers, values)
print(list(res))
Как можно видеть, синтаксис практически тот же, библиотека сама делает всю работу по распараллеливанию вычислений. Впрочем, время выполнения практически то же - 3.4с, ведь мы выполняем вычисления на одном компьютере. При желании читатели могут повторить эксперимент с несколькими компьютерами самостоятельно, описание процесса есть на странице http://ipyparallel.readthedocs.io/en/latest/process.html.
Способ-5
Как говорилось еще в начале книги, Python - это интерпретатор, и код выполняется достаточно медленно. Кардинально повысить скорость можно, используя компиляцию программы перед выполнением, что было реализовано в библиотеке Numba. Ее можно поставить, выполнив команду “pip install numba”.
Метод очень простой и эффективный, и практически не требует изменения кода: всего лишь добавляем перед описанием функции ключевое слово @numba.jit.
Данный код скопирован в отдельный файл Appendix3_test2.py:
import numba
@numba.jit
def get_dividers(n):
divs = 0
for i in xrange(1, n+1):
if n % i == 0:
divs += 1
return divs
Запуск ничем не отличается от предыдущего варианта:
import Appendix3_test2 as test2
res = map(test2.get_dividers, values)
print(list(res))
Время выполнения в обычном, однопоточном режиме - 1.5с. Это уже заметная разница!
И наконец, повторяем вышеописанные процедуры, но функцию заменяем на новый вариант.
Многопоточность:
p = multiprocessing.Pool(processes=8)
res = p.map(test2.get_dividers, values)