Python modul PyCUDA: funkcije i tipovi podatka dostupni u CUDA bibliotekama

Kompletna dokumentacija svih funkcija koje opisujemo u nastavu dostupna je u PyCUDA dokumentaciji, u dijelu koji opisuje funkcije za rad s višedimenzionalnim poljima na GPU-u.

Modul pycuda.gpuarray

Pored dosad korištenih modula, koristimo i pycuda.gpuarray; on radi konverziju iz numpy polja u polja s kojima GPU može manipulirati. Konverzija se inače izvodi automatski korištenjem pomoćnih funkcija pycuda.driver.In(), pycuda.driver.Out() i pycuda.driver.InOut().

Inicijalizacija, izračun rezultata na GPU-u za zbroj matrica je oblika

a = np.ones((20, 20), dtype=np.float32)
b = np.ones((20, 20), dtype=np.float32)

a_gpu = ga.to_gpu(a_cpu)
b_gpu = ga.to_gpu(b_cpu)

result_gpu = ga.empty_like(a_gpu)

zbroji_matrice(result_gpu, a_gpu, b_gpu, block=(20,20,1), grid=(1,1))

Uočite da su result_gpu, a_gpu i b_gpu već gpuarray polja i nema potrebe za pycuda.driver.In(), pycuda.driver.Out() i pycuda.driver.InOut().

Korištenjem činjenice da gpuarray preopterećuje operator zbrajanja, gornji kod može se i dodatno pojednostaviti.

a = np.ones((20, 20), dtype=np.float32)
b = np.ones((20, 20), dtype=np.float32)

a_gpu = ga.to_gpu(a_cpu)
b_gpu = ga.to_gpu(b_cpu)

result_gpu = a_gpu + b_gpu

Modul pycuda.cumath: funkcije koje rade element-po-element na poljima tipa GPUArray

CUDA Math API ima predefinirane funkcije uređaja s jednostrukom preciznošću i dvostrukom preciznošću. Kada se koristi PyCUDA, navedene funkcije su dostupne putem modula pycuda.cumath koji opisujemo u nastavku.

Zaokruživanje i apsolutna vrijednost:

pycuda.cumath.fabs(array)
pycuda.cumath.ceil(array)
pycuda.cumath.floor(array)

Exponencijalne i logaritamske funkcije, korjenovanje:

pycuda.cumath.exp(array)
pycuda.cumath.log(array)
pycuda.cumath.log10(array)
pycuda.cumath.sqrt(array)

Trigonometrijske funkcije:

pycuda.cumath.sin(array)
pycuda.cumath.cos(array)
pycuda.cumath.tan(array)
pycuda.cumath.asin(array)
pycuda.cumath.acos(array)
pycuda.cumath.atan(array)

Računamo li zbroj oblika \(sin(a) + cos(b)\) element po element, kod je oblika

a_cpu = np.ones((20, 20), dtype=np.float32)
b_cpu = np.ones((20, 20), dtype=np.float32)

result_cpu = a_cpu + b_cpu

print("CPU rezultat\n", result_cpu)

a_gpu = ga.to_gpu(a_cpu)
b_gpu = ga.to_gpu(b_cpu)

result_gpu = pycuda.cumath.sin(a_gpu) + pycuda.cumath.cos(b_gpu)

print("GPU rezultat\n", result_gpu.get())

Zadatak

Promijenite kod da računa element po element rezultat oblika \(\sqrt{e^a + \tan(b)}\).

Zadatak

Incijalizirate GPUArray sa 10 redaka i 5 stupaca u kojem ćemo retke tretirati kao vektore. Napišite kod koji vrši zbrajanje svakog vektora sa svakim vektorom a rezultat vraća kao GPUArray veličine \(\binom{10}{2}\) redaka i 5 stupaca.

Aproksimativno računanje broja \(\pi\) kao sume niza

Broj \(\pi\) može se aproksimirati formulom

\[ \pi = \int_0^1 \frac{4}{1 + x^2} dx \approx \frac{1}{n} \sum_{i = 0}^{n - 1} \frac{4}{1 + (\frac{i + 0.5}{n})^2}. \]

sekvencijalni kod je oblika

import numpy as np

n_iterations = 10

def compute_pi(n):
    h = 1.0 / n
    s = 0.0
    for i in range(n):
        x = h * (i + 0.5)
        s += 4.0 / (1.0 + x**2)
    return s * h

pi = compute_pi(n_iterations)
error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

vektorizirani kod koji koristi numpy je oblika

import numpy as np

n_iterations = 1000

def compute_pi(n):
    h = 1.0 / n
    s = 0.0 # nepotrebno
    x = h * (np.arange(1, n) + 0.5)
    s = np.sum(4.0 / (1.0 + x**2))
    return s * h

pi = compute_pi(n_iterations)
error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

paralelni kod je oblika

import pycuda.autoinit
import pycuda.gpuarray as ga
import numpy as np

n_iterations = 1000

def compute_pi(n):
    h = 1.0 / n
    s = 0.0 # nepotrebno
    x = h * (ga.arange(1, n, dtype=np.float32) + 0.5)
    s = ga.sum(4.0 / (1.0 + x**2), dtype=np.float32)
    return s.get() * h

pi = compute_pi(n_iterations)
error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

Zadatak

Modificirajte gornji primjer tako da dodate kod koji mjeri vrijeme izvođenja.
Usporedite vrijeme izvođenja algoritma na CPU-u i na GPU-u za 10⁴, 10⁵, 10⁶ iteracija; ukoliko neku od ovih veličina ne budete mogli izvesti zbog nedostatka memorije na GPU-u, zanemarite tu i sve veće. Opišite svoje zaključke.

Modul pycuda.curandom: generiranje polja sa slučajnim brojevima

pycuda.curandom.rand(shape, dtype=numpy.float32) -- generira GPUArray oblika shape sa "nekvalitetnim" slučajnim brojevima (njihova slučajnost se može dovesti u pitanje)

nekvalitetni_slucajni = pycuda.curandom.rand(10000)
print(nekvalitetni_slucajni)

XORWOW i MRG32k3a generatori

kvalitetniji generatori slučajnih brojeva; oba generatora imaju period barem 2¹⁹⁰
pycuda.curandom.XORWOWRandomNumberGenerator() -- konstruira instancu XORWOW generatora
pycuda.curandom.MRG32k3aRandomNumberGenerator() -- konstruira instancu MRG32k3a generatora

identično sučelje za pristup funkcionalnosti, primjer sa XORWOW generatorom

kvalitetni_slucajni = pycuda.gpuarray.GPUArray(10000, dtype=np.float64)

rng = pycuda.curandom.XORWOWRandomNumberGenerator()

rng.fill_uniform(kvalitetni_slucajni)
print("Uniformna distribucija", kvalitetni_slucajni)

rng.fill_normal(kvalitetni_slucajni)
print("Normalna distribucija", kvalitetni_slucajni)

Zadatak

Promijenite gornji primjer da koristi MRG32k3a generator umjesto XORWOW generatora i da se generiraju 32-bitni brojevi s pomičnim zarezom umjesto 64-bitnih.
Usporedite u terminu vremena izvođenja numpy generator koji se izvodi na CPU-u (funkcija np.random.random()) i PyCUDA MRG32k3a generator koji se izvodi na GPU-u kod generiranja:
- 10000 slučajnih brojeva,
- 100000 slučajnih brojeva,
- 1000000 slučajnih brojeva,
- 10000000 slučajnih brojeva,
- 100000000 slučajnih brojeva.

Tip

Poziv funkcije set_printoptions() iz modula numpy oblika np.set_printoptions(threshold=np.nan) čini da se ispisuju svi brojevi numpy polja, bez obzira na njegovu veličinu.

Generatori zasnovani sa Soboljevim nizovima

Info

Generatori Sobol32, ScrambledSobol32, Sobol64 i ScrambledSobol64 zasnovani su na Soboljevim nizovima. Više informacija o njima možete naći u službenoj PyCUDA dokumentaciji.

Aproksimativno računanje broja \(\pi\) korištenjem Monte Carlo metode

Monte Carlo metode su skupina metoda koja na temelju slučajnih brojeva i velikog broja pokusa daju aproksimaciju određene vrijednosti.

Broj \(\pi\) može se aproksimativno računati i korištenjem Monte Carlo metode. Da bi to vidjeli, uzmimo da je u koordinatnom sustavu ucrtan jedinični krug unutar jediničnog kvadrata. Promotrimo njegovu četvrtinu u prvom kvadrantu:

četvrtina kvadrata je površine 1 (obzirom da su mu obje stranice duljine 1),
čevrtina kruga je površine \(\frac{\pi}{4}\) (obzirom da mu je radijus 1).

Slučajno odabrana točka unutar četvrtine kvadrata ima vjerojatnost \(\frac{\pi}{4}\) da upadne unutar četvrtine kruga. Dakle, ako s \(n\) označimo broj slučajno odabranih točaka, a s \(h\) broj točaka koje se od \(n\) slučajno odabranih nalaze unutar četvrtine kruga, aproksimativno možemo odrediti broj \(\pi\) kao

\[ \pi = 4 \times \frac{h}{n}. \]

Povećavajući \(n\) dobivamo točniju aproksimaciju, a navedena metoda naziva se Monte Carlo metoda. Programski kod koji aproksimira broj \(\pi\) korištenjem Monte Carlo metode je oblika:

sekvencijalni kod

import numpy as np

n_random_choices = 10000
hits = 0
throws = 0

for i in range (0, n_random_choices):
    throws += 1
    x = np.random.random()
    y = np.random.random()
    dist = np.math.sqrt(x * x + y * y)
    if dist <= 1.0:
        hits += 1

pi = 4 * (hits / throws)

error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

paralelni kod

lošiji pristup; for petlja vrši 10000 * 4 dohvaćanja elemenata iz memorije GPU-a

import pycuda.autoinit
import pycuda.driver as drv
import pycuda.gpuarray as ga
import pycuda.curandom
import pycuda.cumath
from pycuda.compiler import SourceModule
import numpy as np

n_random_choices = 10000
hits = 0
throws = n_random_choices

x_polje = pycuda.curandom.rand(n_random_choices)
y_polje = pycuda.curandom.rand(n_random_choices)
for i in range (0, n_random_choices):
    dist = np.math.sqrt(x[i] * x[i] + y[i] * y[i])
    if dist <= 1.0:
        hits += 1

pi = 4 * (hits / throws)

error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

bolji pristup, eliminacija for petlje, radimo na polju

import pycuda.autoinit
import pycuda.driver as drv
import pycuda.gpuarray as ga
import pycuda.curandom
import pycuda.cumath
from pycuda.compiler import SourceModule
import numpy as np

n_random_choices = 10000
hits = 0
throws = n_random_choices

x_polje = pycuda.curandom.rand(n_random_choices)
y_polje = pycuda.curandom.rand(n_random_choices)
dist = pycuda.cumath.floor(pycuda.cumath.sqrt(x_polje * x_polje + y_polje * y_polje))
hits = n_random_choices - dist.get().sum()

pi = 4 * (hits / throws)

error = abs(pi - np.math.pi)
print("pi is approximately %.16f, error is approximately %.16f" % (pi, error))

Zadatak

Modificirajte primjer dan za Monte Carlo simulaciju da dodate kod koji mjeri vrijeme izvođenja.
Usporedite vrijeme izvođenja simulacije na CPU-u i na GPU-u za 10⁴, 10⁵, 10⁶, 10⁷, 10⁸, 10⁹ iteracija; ukoliko neku od ovih veličina ne budete mogli izvesti zbog nedostatka memorije na GPU-u, zanemarite tu i sve veće. Opišite svoje zaključke.

Zadatak

Usporedite vrijeme izvođenja programa i točnost aproksimacije pi korištenjem:

sumiranja elemenata niza,
Monte Carlo simulacije,

na CPU-u i na GPU-u za 10⁴, 10⁵, 10⁶ iteracija. Opišite svoje zaključke.

Zadatak

Simulirajmo ekonomiju. Inicijalizirajte GPUArray veličine 200 elemenata koji predstavlja financijsko stanje 200 jediniki u ekonomiji koju simuliramo. Incijalizirajte vrijednosti elemenata na 10. Definirajte dvije operacije

porast_pad_ekonomije(), koji svačije financijsko stanje množi sa koeficijentom k koji se određuje kao slučajan broj u intervalu \([0.5, 1.5]\).
preraspodjela_trgovina(), koja vrši transakciju određenog iznosa (slučajan broj u rasponu \([0, 0.5]\)) između dvije slučajno odabrane jedinke.

Izvedite Monte Carlo simulaciju veličine 100000 iteracija od kojih svaka iteracija vrši

slučajno između 5 i 10 operacija preraspodjela_trgovina
jednu operaciju porast_pad_ekonomije.

Author: Vedran Miletić