5. Если простое число найдено, мы кладем его в глобальный массив results. Но т.к. несколько ядер работают с одним массивом одновременно, мы должны быть уверены что данные будут сохранены корректно, для этого мы используем функцию atomic_add. Первую ячейку массива мы таким образом используем для хранения количества найденных простых чисел.
6. Когда выполнение программы завершено, мы загружаем данные обратно с помощью функции dest_dev.get(). Т.к. ядра запускались параллельно, то простые числа лежат в массиве вперемешку, чтобы их отсортировать, мы вызываем функцию numpy.sort.
Как можно видеть, процесс довольно-таки громоздкий, но оно того стоит. Однопоточная программа, написанная на Си с таким же кодом, выполнялась 14 секунд. Вышеприведенная программа, запускаемая на видеокарте, выполнила расчет за 1.2 секунды.
Разумеется, еще раз стоит повторить, что “игра стоит свеч” лишь в том случае, если задача хорошо распараллеливается на небольшие блоки, в таком случае выигрыш будет заметен.
Рассмотрим теперь тот же код, но написанный с помощью библиотеки pycuda. Разумеется, выполнить его смогут только владельцы видеокарт NVIDIA.
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
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и.