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:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
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")
In [3]:
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")
In [4]:
(free,total)=cuda.mem_get_info()
free,total
Out[4]:
Ponieważ chcemy zapisywać liczny w pojedynczej prezycji, to mamy możliwość zaalokowania pamięci na:
In [5]:
free/4
Out[5]:
liczb 4 bajtowych.
In [6]:
%%time
print(cuda.mem_get_info()[0]/1024**2)
N = 500000000
a = curand((N//2,))
b = curand((N//2,))
ctx.synchronize()
print(cuda.mem_get_info()[0]/1024**2)
In [7]:
a.get().shape
Out[7]:
In [8]:
%%time
for i in range(10):
a1 = a.get()
In [9]:
cuda.mem_get_info()[0]/1024.**3
Out[9]:
In [10]:
%%time
pole = float(krnl(a, b).get())/(N/2)
ctx.synchronize()
print( "OK")
In [11]:
print (4*pole,"Pi z %d milionów losowań!"%(N/2/1e6))
print (np.pi)
In [ ]:
In [ ]:
In [12]:
x,y = a.get()[::10000],b.get()[::10000]
plt.plot(x,y,'.')
c = (x-0.5)**2+(y-0.5)**2<0.25
plt.plot(x[c],y[c],'.')
Out[12]:
In [13]:
try:
ctx.pop()
ctx.detach()
print ("OK!")
except:
print( "No CTX!")
In [14]:
%load_ext Cython
In [15]:
%%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(5000000):
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
In [16]:
%%time
cpu_rand()
Out[16]:
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$:
In [17]:
628317563/800000000.*4
Out[17]:
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.
In [ ]: