Skirtumai

Čia matote skirtumus tarp pasirinktos versijos ir esamo dokumento.

--- hpc:gpu [2022/01/31 16:46] – linp
+++ hpc:gpu [2022/01/31 17:11] – linp
@@ Linija 77: / Linija 77: @@
 </code>
+====== GPU išnaudojimas lygiagretinant kodą su Numba biblioteka ======
+=== Įvadas į CUDA Python su Numba  ===
+CUDA skaičiavimo biblioteka įgalina programų pagreitinimą, vykdant skaičiavimus ant GPU.
+**Numba** yra "vykdymo-realiu-laiku" Python funkcijų kompiliatorius, Python funkcijoms paspartinti. **Numba** dažnai naudojama Python programuotojų norint, norintiems GPU paspartinti savo programas, jų neperašant į C/C++ kodą, ypač jei kodas jau naudoja **NumPy** masyvuos.
+=== Pavyzdžiai darbas su Numba: Pvz. 1  ===
+Pavydžiui, realizuokime Pitagorinę sudėtį:
+https://en.wikipedia.org/wiki/Hypot
+<code shell>
+# importuojame jit kompiliatorių
+from numba import jit
+import numpy as np
+import math
+# Sintaksė @jit ekvivalentu užrašui `hypot = jit(hypot)`.
+@jit
+def hypot(x, y):
+    x = abs(x);
+    y = abs(y);
+    t = min(x, y);
+    x = max(x, y);
+    t = t / x;
+    return x * math.sqrt(1+t*t)
+# Bandome
+hypot(3.0, 4.0)
+# Nesukompiliuotos su *jit* funkcijos iškvietimas
+ypot.py_func(3.0, 4.0)
+</code>
+Gauti metodų vykdymo laikus galėtume atitinkamai pvz:
+<code shell>
+timeit hypot.py_func(3.0, 4.0)
+timeit hypot(3.0, 4.0)
+</code>
+=== Pavyzdžiai darbas su Numba: Pvz. 2  ===
+Parašykime programą pi vertinimui Monte Karlo metodu (žr. https://academo.org/demos/estimating-pi-monte-carlo/)
+<code shell>
+from numba import jit
+import random
+@jit
+def monte_carlo_pi(nsamples):
+    acc = 0
+    for i in range(nsamples):
+        x = random.random()
+        y = random.random()
+        if (x**2 + y**2) < 1.0:
+            acc += 1
+    return 4.0 * acc / nsamples
+import matplotlib.pyplot as plt
+nsamples = 1000000
+n = [2**i for i in range(0, 20)]
+pi_values = [monte_carlo_pi(i) for i in n]
+plt.plot(n, pi_values)
+plt.axhline(y=np.pi, color='r', linestyle='-')
+plt.xscale('log')
+plt.xlabel("Bandymų skaičius n")
+plt.ylabel("pi įvertis")
+</code>
+Gauti metodų vykdymo laikus galėtume atitinkamai pvz:
+<code shell>
+timeit monte_carlo_pi(nsamples)
+timeit monte_carlo_pi.py_func(nsamples)
+</code>
+=== Kaip Numba veikia?  ===
+Numba kompiliatorius atsivželgia į duomenų tipus ir optimizuoja tarpinius skaičiavimu (angl. Intermediate representation (IR))
+### **Numba** optimizuojant kodą skirtą GPU naudojant **NumPy** universalias funkcijas (ufuncs)
+Toliau laikysime GPU programavimo aprėpyje jog dirbsime su NumPy universaliomis funkcijas (arba ufuncs) - t.y. standartinės bibliotekos funkcijos žr. dokumentaciją:   https://docs.scipy.org/doc/numpy-1.15.1/reference/ufuncs.html
+=== Pavyzdžiai darbas su Numba: Pvz. 3  ===
+<code shell>
+import numpy as np
+a = np.array([1, 2, 3, 4])
+b = np.array([10, 20, 30, 40])
+np.add(a, b)
+# Universalios funkcijos veikia su vektoriais ir skaičiais
+np.add(a, 100)
+# Matriciniu atveju, vektoriai taikomi sąrašo elementams t.y. eilutėms
+c = np.arange(4*4).reshape((4,4))
+print('c:', c)
+np.add(b, c)
+</code>
+=== Universalių funkcijų kūrimas ir optimizavimas  ===
+Svarbu suprasti, jog norint kurti universalias fukcijas, tokias turime ir naudoti. Šiuo atveju, tam mum bus reikalinkas metodas *vectorize*.
+Šiame pačiame pirmame pavyzdyje naudosime *vectorize*, kad sudarytume ir optimizuotume sukurtą universalią funkciją **CPU**
+<code shell>
+from numba import vectorize
+@vectorize
+def add_ten(num):
+    return num + 10
+nums = np.arange(10)
+nums
+add_ten(nums)
+# Dabar norint išnaudoti **GPU** optimizavimą mum reikia
+@vectorize(['int64(int64, int64)'], target='cuda')
+def add_ufunc(x, y):
+    return x + y
+add_ufunc(a, b)
+</code>
+Gauti metodų vykdymo laikus galėtume atitinkamai pvz:
+<code shell>
+timeit np.add(b, c)   # NumPy - CPU
+timeit add_ufunc(b, c) # Numba - GPU
+</code>
+# Ant CPU veikia greičiau :)
+# 1. duomenys per maži
+# 2. operacijos primityvios
+# 3. kopijuojame duomenis į GPU
+# 4. naudojame didelius duomenis (int64)
+=== Pavyzdžiai darbas su Numba: Pvz. 4  ===
+Panagrinėkime sudėtingesnį atvejį f(x|\mu, \sigma) = \frac{e^{-((x-\mu)/\sigma)^2}}{\sigma \sqrt{2\pi}}
+<code shell>
+import math
+SQRT_2PI = np.float32((2*math.pi)**0.5)
+@vectorize(['float32(float32, float32, float32)'], target='cuda')
+def gaussian_pdf(x, mean, sigma):
+    return math.exp(-0.5 * ((x - mean) / sigma)**2) / (sigma * SQRT_2PI)
+import numpy as np
+x = np.random.uniform(-3, 3, size=1000000).astype(np.float32)
+x = np.sort(x)
+mean = np.float32(0.0)
+sigma = np.float32(1.0)
+f = gaussian_pdf(x, 0.0, 1.0)
+plt.plot(x, f)
+plt.xlabel('x')
+plt.ylabel('f(x)')
+</code>
+Gauti metodų vykdymo laikus galėtume atitinkamai pvz:
+<code shell>
+timeit gaussian_pdf(x, mean, sigma)
+import scipy.stats
+norm_pdf = scipy.stats.norm
+timeit norm_pdf.pdf(x, loc=mean, scale=sigma)
+</code>
+=== Pavyzdžiai darbas su Numba: Pvz. 5  ===
+Tais atvejais, kai norime spartinti vektorinius skaičiavimus, o ne skaičiavimus paelemenčiui -- naudoti *numba.cuda.jit*
+<code shell>
+from numba import cuda
+@cuda.jit(device=True)
+def polar_to_cartesian(rho, theta):
+    x = rho * math.cos(theta)
+    y = rho * math.sin(theta)
+    return x, y
+@vectorize(['float32(float32, float32, float32, float32)'], target='cuda')
+def polar_distance(rho1, theta1, rho2, theta2):
+    x1, y1 = polar_to_cartesian(rho1, theta1)
+    x2, y2 = polar_to_cartesian(rho2, theta2)
+    return ((x1 - x2)**2 + (y1 - y2)**2)**0.5
+n = 1000000
+rho1 = np.random.uniform(0.5, 1.5, size=n).astype(np.float32)
+theta1 = np.random.uniform(-np.pi, np.pi, size=n).astype(np.float32)
+rho2 = np.random.uniform(0.5, 1.5, size=n).astype(np.float32)
+theta2 = np.random.uniform(-np.pi, np.pi, size=n).astype(np.float32)
+polar_distance(rho1, theta1, rho2, theta2)
+</code>
+=== Plačiau apie GPU išnaudojimą  ===
+ * [[https://klevas.mif.vu.lt/~linp/hpc/ | Medžiaga ir pavyzdžiai: klevas.mif.vu.lt/~linp/hpc/ ]]
+ * [[https://drive.mif.vu.lt/s/9Zpzxzq4XCyMozp|GPU skaičiavimai HPC insfrastruktūroje (video) ]]