Pierwszy program wykorzystujący CUDA, będzie obliczał $\pi$ metodą statystyczną.
Jedną z metod obliczania przybliżonej wartości liczby $\pi$ jest metoda polegająca na pomiarze prawdopodobieństwa trafiena o obszar, będący kołem. Załóżmy, że mamy możliwość wylosowania z rozkładem jednostajnym punktu w kwadracie o jednostkowym boku. Jeżeli teraz zliczymy ilośc punktów zawartych w dowolnej figurze wewnątrz tego kwadratu, to liczba ta będzie proporcjonalna do powierzfchni tej figury. Zastosujmy tą obserwację do koła o promieniu $r = \frac{1}{2}$ wpisanego w nasz kwadrat. Powierznia takiego koła wynosi $P=\pi r^2 = \frac{\pi}{4}$. Obliczając więc stosunek liczby punktów zawartych w tym kole do wszystkich punktów otrzymamy $\frac{\pi}{4}$.
Wykorzystamy bibliotekę pyCUDA w w szczególnosci funkcjonalność gpuarray. Ponadto będziemy korzystać z wbudowanych generatorów liczb losowych.
Wykorzystanie gpuarray nie wymaga zbyt obszernej wiedzy o architekturze CUDA. Nie będziemy musieli pisać własnego "kernela". Pomino tego gpuaray oferuje niezwykle dużo możliwości, które mogą być przydatne do rozwiązywania problemów napotkanych w naukach ścisłych lub inżynierii.
Wysumowanie wartości w wektorze wykonane z wykorzystaniem procesora równoleglego jest zwane problemem redukcji równoległej. Efektywna implementacja redukcji równoległej na CUDA wymaga mistrzowskiej znajomości architektury. Naiwne implementacje z reguły zaniżająpotencjalną wydajność kilka lub kilkadziesiąt razy. W bibliotece pycuda mamy zaimplementowany szablo jądra redukcji, który możemy wykorzytać, nie wnikając w szczegóły jego implementacji.
Nasz program będzie składał się z następujących elementów:
Należe podkreślić, że taka metoda obliczania liczby $\pi$ jest wyjątkowo kiepska. Dlateczogo więc chcemy ja zaimplementować i to do tego na GPU?
Poniższa implementacja ma następujące walory dydaktyczne:
import pycuda.gpuarray as gpuarray
import numpy
from pycuda.curandom import rand as curand
from pycuda.compiler import SourceModule
import pycuda.driver as cuda
try:
ctx.pop()
ctx.detach()
except:
print "No CTX!"
cuda.init()
device = cuda.Device(0)
ctx = device.make_context()
print device.name(), device.compute_capability(),device.total_memory()/1024.**3,"GB"
print "W systemie mamy :",cuda.Device.count(), " urządzenia"
GeForce GTX 680 (3, 0) 3.99932861328 GB W systemie mamy : 1 urządzenia
Kod operacji wykonanej na każdej parze liczb przed redukcją:
signbit( (powf(x[i]-0.5f,2.0f)+powf(y[i]-0.5f,2.0f))-0.25f )
Zauważmy:
from pycuda.reduction import ReductionKernel
krnl = ReductionKernel(numpy.dtype(numpy.int32), neutral="0",
reduce_expr="a+b",
map_expr="signbit( (powf(x[i]-0.5f,2.0f)+powf(y[i]-0.5f,2.0f))-0.25f )",
arguments="float *x, float *y")
Funkcja curand pozwala nam skorzystać z wbudowanego generatora liczb losowych w CUDA. Możemy z jej pomocą zapełnić całą pamięć liczbami losowymi z przedziału $(0,1)$. Sprawdźmy najpiew ile mamy na GPU pamięci:
(free,total)=cuda.mem_get_info()
free,total
(3515383808, 4294246400)
Ponieważ chcemy zapisywać liczny w pojedynczej prezycji, to mamy możliwość zaalokowania pamięci na:
free/4
878845952
liczb 4 bajtowych.
%%time
print cuda.mem_get_info()[0]/1024**2
N = 800000000
a = curand((N/2,))
b = curand((N/2,))
ctx.synchronize()
print cuda.mem_get_info()[0]/1024**2
3352 300 CPU times: user 196 ms, sys: 56 ms, total: 252 ms Wall time: 420 ms
Poeksperymentuj ze zwalnianiem pamięci wykonując np. takie operacje:
del a
print cuda.mem_get_info()[0]/1024**2
del b
print cuda.mem_get_info()[0]/1024**2
%%time
for i in range(10):
a1= a.get()
CPU times: user 3.62 s, sys: 3.84 s, total: 7.46 s Wall time: 7.48 s
cuda.mem_get_info()[0]/1024.**3
0.290557861328125
%%time
pole = float(krnl(a, b).get())/(N/2)
ctx.synchronize()
print "OK"
OK CPU times: user 96 ms, sys: 0 ns, total: 96 ms Wall time: 95.3 ms
print 4*pole,"Pi z %d milionów losowań!"%(N/2/1e6)
print np.pi
3.14158223 Pi z 400 milionów losowań! 3.14159265359
try:
ctx.pop()
ctx.detach()
print "OK!"
except:
print "No CTX!"
OK!
... warto uważać ile punktów się rysuje...
x,y = a.get()[::100000],b.get()[::100000]
plot(x,y,'.')
c = (x-0.5)**2+(y-0.5)**2<0.25
plot(x[c],y[c],'.')
[<matplotlib.lines.Line2D at 0x364ead0>]
Aby się przekonać, czy GPU rzeczywiście wykonał powyższe operacje szybciej niż CPU, sprawdźmy jak szybko wykona się kod skompilowany przez cythona. Cython daje wyniki nie gorsze od czystego "C".
%load_ext cythonmagic
The cythonmagic extension is already loaded. To reload it, use: %reload_ext cythonmagic
%%cython
cdef extern from "stdlib.h":
int RAND_MAX
from libc.stdlib cimport rand
from libc.math cimport pow
def cpu_rand():
cdef double a,b
cdef long n=0
for i in range(800000000):
a = float(rand())/RAND_MAX
b = float(rand())/RAND_MAX
if pow(a-0.5,2.0)+pow(b-0.5,2.0)-0.25<0:
n=n+1
return n
%%time
cpu_rand()
CPU times: user 26.7 s, sys: 0 ns, total: 26.7 s Wall time: 26.6 s
628317563
Mamy więc ok 30s na CPU i ok 300ms na GPU - czyli ok 100x ! Należy jednak zauważyć, że na GPU wykonywaliśmy alokację wszystkich liczb w pamięci a potem ich równoległą redukcję. Jest to niezbyt optymalne i można by się spodziewać dalszego wzrostu względnej wydajności jesli by porównywać dokładnie te same algorytmy.
Zobaczmy jeszcze czy na CPU $\pi=3.14$:
628317563/800000000.*4
3.141587815
Czy można zmodyfikować algorytm tak by obiczał $\pi$ jeszcze szybciej? (oczywiście tą samą metodą). Może warto by z pomocą jednego wątku generować nie jedną liczbę losową ale np 10, lub 100 i obliczać sumę częściową?
Wykonać identyczny funkcjonalnie algorytm na CPU. Użyć C/C++ lub cythona.