import numpy, time
# Define kernel
mod = SourceModule('''
__global__ void primes(int *results, int maxVal) {
int gid = blockDim.x*blockIdx.x + threadIdx.x;
if (gid < 2 || gid >= maxVal) return;
int num = gid + 1;
for(int p=2; p<=num/2 + 1; p++) {
if ((num % p) == 0) {
return;
}
}
int val = atomicAdd(&results[0], 1);
results[val+1] = num;
}''')
N = 500000
n, grid = 1 + N/512, 512
# Host array
h_array = numpy.zeros(n*grid, dtype=numpy.int32)
# Copy host array to device array
d_array = cuda.mem_alloc(h_array.nbytes)
cuda.memcpy_htod(d_array, h_array)
# Invoke kernel
func = mod.get_function("primes")
func(d_array, numpy.int32(N), block=(n, 1, 1), grid=(grid, 1))
# Copy device array to host array
cuda.memcpy_dtoh(h_array, d_array)
count = h_array[0]
primes = numpy.sort(h_array[1:count+1])
print "Found prime numbers:", count
# for p in primes:
# print p
Как можно видеть, принцип тот же. Единственная разница состоит в том, что “вручную” пришлось указать размер блока в 512 значений, и весь массив разбить на блоки, используя параметр grid - библиотека OpenCL делала этот процесс автоматически. И разумеется, синтаксис и названия функций разных библиотек тоже разные, хотя суть остается примерно та же.
Приложение 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
Быстрее, но лишь чуть-чуть.
Заметно лучший результат можно получить, используя профессиональную математическую библиотеку numpy, хотя она довольно сложна для новичков.
import numpy as np
start_time = time.time()
a = np.arange(1,1000001, dtype=np.uint64)
s = np.sum(a ** 2)
print("Sum = {}, T = {}s".format(s, time.time() - start_time))
Результат 0.008с - в 40 раз быстрее “обычного” Python-кода.
И наконец, призываем “тяжелую артиллерию”: перепишем программу на языке 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-библиотек - для многих задач такие библиотеки уже написаны на языке C++, и их можно использовать, и они работают довольно-таки быстро.
Для тех, кто захочет максимально широко использовать Python в математических расчетах, крайне рекомендуется ознакомиться с библиотекой numpy, она является стандартом де-факто для всех серьезных расчетов.
Теперь посмотрим, как еще можно ускорить Python-программу с помощью многопоточности и других хитростей.
Приложение 3 - Пример ускорения Python-программы
Рассмотрим практический пример сравнения скорости вычислений на 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, позволяющий разбить вычисления над массивом на несколько потоков (способ работает только в Python 3).
if __name__ == '__main__':
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
Запуск ничем не отличается от предыдущего варианта:
res = map(test2.get_dividers, values)
print(list(res))
Время выполнения в обычном, однопоточном режиме - 1.5с. Это уже заметная разница!
И наконец, повторяем вышеописанные процедуры, но функцию заменяем на новый вариант.
Многопоточность:
p = multiprocessing.Pool(processes=8)
res = p.map(test2.get_dividers, values)
Время выполнения: 1.04c.
И последний вариант, библиотека ipyparallel:
res = dv.map_sync(test2.get_dividers, values)
print(list(res))
Время выполнения: 0.57c.
Таким образом, по сравнению с изначальным вариантом, мы увеличили скорость расчетов почти в 20 раз, при этом сам код изменять или переписывать практически не пришлось.
В целом, можно предложить следующий алгоритм оптимизации вычислений на языке Python:
1) По возможности, стараться использовать встроенные функции, имеющиеся в библиотеках math, numpy и др. Они уже написаны на C или C++, и работают максимально эффективно.
2) Использовать библиотеку numba, указав в начале функции префикс @numba.jit. Это самый простой и “безболезненный” способ, он чуть увеличит время запуска программы, но скорость расчетов может вырасти в 10 раз. Правда, этот способ не поможет, если код уже состоит только из вызовов встроенных функций из п1. Если же код содержит циклы, условия, другие проверки, то скорость выполнения вырастет значительно.
3) Если скорость расчета нужно еще увеличить, переходим к более тяжелой артиллерии - многопоточности. Используя библиотеку multiprocessing, можно разбить вычисления на несколько ядер процессора - это может увеличить скорость в 2-4 раза. Но важно иметь в виду, что межпроцессное взаимодействие - сама по себе довольно-таки долгая операция, так что делать ее имеет смысл для действительно долго выполняющихся фрагментов кода. Небольшая функция, вызывающаяся много раз в разных потоках, будет работать в итоге даже медленнее за счет дополнительных “накладных расходов”.
4) И наконец, если скорость расчета нужно еще увеличить, можно разбить вычисления на несколько компьютеров в сети, используя библиотеку ipyparallel. Как и в предыдущем случае, заранее следует продумать структуру алгоритма, чтобы производить вычисления крупными блоками, это сделает метод более эффективным.
5) Если необходимо обработать большой объем данных, целесообразно рассмотреть вычисления на видеокарте (GPU). За счет параллельности обработки, можно запустить программу одновременно на нескольких сотнях ядер видеокарты, что может дать огромный прирост в скорости. Такие библиотеки (OpenCL и CUDA) есть практически под все видеокарты, имеющиеся сегодня в продаже.