Norint pasinaudoti PST skaičiavimo resursais su GPU reikės formuoti užduočių skriptus (.sh):

Pavyzdys, bandant patikrinti ar jum suteikta GPU. Sukuriamas failas gpu.txt su indikatoriumi ir vykdymo laiko momentu:

test_gpu.sh

#!/bin/sh
#SBATCH -p gpu
#SBATCH -n1
#SBATCH --gres gpu
 
source gpu_env/bin/activate
python3 test.py

Pasinagrinėkime pavyzdį paeilučiui:

$ #SBATCH -p gpu 

nukreipiama į PST su vaizdo plokščių resursais skaičiavimo resursą.

$ #SBATCH -n1

nurodome, koks CPU poreikis bus reikalingas

$ #SBATCH --gres gpu

nurodome, koks GPU poreikis bus reikalingas (esant N GPU poreikiui būtų nurodoma gpu:N).

$ source gpu_env/bin/activate

pasiruošiame aplinką darbui (pirmą kartą ją reikia sukurti atskirai žr. aplinka.sh)

$ python3 test.py

Kviečiame savo kodo skriptą, atlikti užduočiai su Python

Čia iškviečiamo Python kalbos skripto kodas, iškviečiant Pytorch biblioteką ir patikrinat GPU prieigą, kurios dažniausiai reikia kuriant mašininio mokymo modelius.

test.py

import torch
import time
 
val = torch.cuda.is_available()
 
f = open("gpu.txt", "a")
ts = time.time()
f.write("GPU is loaded {} at {}\n".format(val, ts))
f.close()

Čia iškviečiamo Python kalbos skripto kodas, iškviečiant Pytorch biblioteką ir patikrinat GPU prieigą, kurios dažniausiai reikia kuriant mašininio mokymo modelius.

aplinka.sh

#!/bin/bash
 
python3 -m venv gpu_env
 
# pirmą kartą susirašome reikiamas bibliotekas savo projektui pvz:
 
pip3 install wheel
pip3 install pillow scikit-image
pip3 install numba 
pip3 install torch==1.10.0+cu113 torchvision==0.11.1+cu113 torchaudio==0.10.0+cu113 -f https://download.pytorch.org/whl/cu113/torch_stable.html

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.

Pavydžiui, realizuokime Pitagorinę sudėtį:

https://en.wikipedia.org/wiki/Hypot

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

Gauti metodų vykdymo laikus galėtume atitinkamai pvz:

timeit hypot.py_func(3.0, 4.0)
 
timeit hypot(3.0, 4.0)

Parašykime programą pi vertinimui Monte Karlo metodu (žr. https://academo.org/demos/estimating-pi-monte-carlo/)

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

Gauti metodų vykdymo laikus galėtume atitinkamai pvz:

timeit monte_carlo_pi(nsamples)
 
timeit monte_carlo_pi.py_func(nsamples)

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

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)

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

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)

Gauti metodų vykdymo laikus galėtume atitinkamai pvz:

timeit np.add(b, c)   # NumPy - CPU
 
timeit add_ufunc(b, c) # Numba - GPU

# Ant CPU veikia greičiau :) # 1. duomenys per maži # 2. operacijos primityvios # 3. kopijuojame duomenis į GPU # 4. naudojame didelius duomenis (int64)

Panagrinėkime sudėtingesnį atvejį f(x|\mu, \sigma) = \frac{e^{-((x-\mu)/\sigma)^2}}{\sigma \sqrt{2\pi}}

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

Gauti metodų vykdymo laikus galėtume atitinkamai pvz:

timeit gaussian_pdf(x, mean, sigma)
 
import scipy.stats 
norm_pdf = scipy.stats.norm
timeit norm_pdf.pdf(x, loc=mean, scale=sigma)

Tais atvejais, kai norime spartinti vektorinius skaičiavimus, o ne skaičiavimus paelemenčiui – naudoti *numba.cuda.jit*

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)

* Medžiaga ir pavyzdžiai: klevas.mif.vu.lt/~linp/hpc/

* GPU skaičiavimai HPC insfrastruktūroje (video)