Obliczamy $\pi$ metodą Monte-Carlo w CUDA.

Pierwszy program wykorzystujący CUDA, będzie obliczał $\pi$ metodą statystyczną.

Opis metody

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}$.

  1. Kiepska metoda obliczania $\pi$
  2. Użycie curand oraz równoległej redukcji w pycuda
  3. Ile czasu zajmuje posumowanie 4GB RAM ;-)

Implementacja

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.

Redukcja równoległa

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.

Struktura programu

Nasz program będzie składał się z następujących elementów:

  1. Zainicjalizujemy kontekst na wybranym urządzeniu
  2. Definiujemy jądro redukcji z szablonu w pyCUDA. Operatorem redukcji będzie suma, a przed wykonaniem operacji wywolujemy funkcję która daje jedynkę jesli punkt należy do naszego koła
  3. Generujemy losowe liczby z przedziału $(0,1)$ bezpośrednio w pamięci urządzenia.
  4. Wykonujemy jądro na liczbach losowych i dzielimy otrzymaną wartość przez $4$. Powinniśmy otrzymać przybliżoną wartość $\pi$.

Co jest ciekawego w tym algorytmie?

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:

  1. Pokazuje naocznie jak szybko przebiegają operacje - np. zapełnienia 1GB pamięci liczbami losowymi
  2. Korzysta z kernela redukcji równoległej. Jest to niezwykle pożyteczny kernel. Warto zauważyć, że wykonuje on de facto operacje identyczną do tej, która występuje w całkowanie numerycznym - czyli liczeniu kwadratury. Innymi słowy mamy szablon programu całkującego!

Inicjalizacja


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")


No CTX!
Tesla K40m (3, 5) 11.1717529296875 GB
W systemie mamy : 2  urządzenia

Jądro redukcji

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:

  1. Używamy funkcji i zmiennych o pojedyńczej precyzji
  2. Funkcja signbit zwraca znak wyrażenia jako liczbę całkowitą

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")

Generacja liczb losowych na urządzeniu

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:


In [4]:
(free,total)=cuda.mem_get_info()
free,total


Out[4]:
(11918901248, 11995578368)

Ponieważ chcemy zapisywać liczny w pojedynczej prezycji, to mamy możliwość zaalokowania pamięci na:


In [5]:
free/4


Out[5]:
2979725312.0

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)


11366.75
9459.25
CPU times: user 108 ms, sys: 28 ms, total: 136 ms
Wall time: 423 ms

Zadanie

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

In [7]:
a.get().shape


Out[7]:
(250000000,)

In [8]:
%%time 
for i in range(10):
    a1 = a.get()


CPU times: user 1.21 s, sys: 1.78 s, total: 2.99 s
Wall time: 2.99 s

In [9]:
cuda.mem_get_info()[0]/1024.**3


Out[9]:
9.237548828125

In [10]:
%%time 

pole = float(krnl(a, b).get())/(N/2)
ctx.synchronize()

print( "OK")


OK
CPU times: user 20 ms, sys: 12 ms, total: 32 ms
Wall time: 31.1 ms

In [11]:
print (4*pole,"Pi z %d milionów losowań!"%(N/2/1e6))
print (np.pi)


3.141682432 Pi z 250 milionów losowań!
3.141592653589793

In [ ]:


In [ ]:

Wykres

... warto uważać ile punktów się rysuje...


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]:
[<matplotlib.lines.Line2D at 0x7f2ac8203c88>]

In [13]:
try:
    ctx.pop()
    ctx.detach()
    print ("OK!")
except:
    print( "No CTX!")


OK!

CPU

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".


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()


CPU times: user 240 ms, sys: 0 ns, total: 240 ms
Wall time: 238 ms
Out[16]:
3928039

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]:
3.141587815

Zadania

Szybsze obliczanie $\pi$ na GPU

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ą?

Algorytm na CPU z alokacją pamięci.

Wykonać identyczny funkcjonalnie algorytm na CPU. Użyć C/C++ lub cythona.


In [ ]: