9. cvičenie – CUDA pomocou Numba

Zdroj: https://nyu-cds.github.io/python-numba/05-cuda/

Úvod do programovania na GPU

Po nainštalovaní modulu numba overte, či váš systém obsahuje grafickú kartu s podporou CUDA. Spustite interpreter jazyka Python a vykonajte tento kód:

Python
from numba import cuda
print(cuda.gpus)

				1
2

						from numba import cuda
print(cuda.gpus)

Na výstupe by ste mali vidieť:

<Managed Device 0>

				1

						<Managed Device 0>

Ak je CUDA funkčná, môžete pokračovať ďalej v cvičení.

Ak sa vyskytla chyba, váš systém nepodporuje CUDA a musíte použiť jej simuláciu. Podľa pokynov v prednáške pred spustením prostredia IDLE alebo iného pracovného prostredia nastavte premennú prostredia NUMBA_ENABLE_CUDASIM na hodnotu 1.

Definícia a spustenie kernel funkcie

Kernel funkcia je typ funkcie, ktorej beh sa síce vyvoláva na CPU, ale jej telo beží na GPU. Charakterizujú ju dve hlavné veci:

Nesmie vracať žiadnu hodnotu. Všetky výsledky výpočtu je nutné odovzdávať pomocou poľa na vstupe funkcie.
Pri spustení kernel funkcie sa definuje štruktúra výpočtových vlákien: počet vlákien v bloku a počet blokov. Takáto parametrizácia umožňuje, aby bol ten istý kernel volaný v rôznych konfiguráciách (podľa možností toho-ktorého zariadenia GPU).

Nasleduje ukážka definície kernelu:

Python
from numba import cuda


@cuda.jit
def my_kernel(io_array):
    """
    Code for kernel.
    """
    # code here

				
					
				1
2
3
4
5
6
7
8
9

						from numba import cuda
 
 
@cuda.jit
def my_kernel(io_array):
    """
    Code for kernel.
    """
    # code here

					

			

Pri spustení je najdôležitejšie vhodne zvoliť parametre výpočtu – počet blokov a počet vlákien v jednom bloku. Celkový počet výpočtových vlákien na GPU je daný ako súčin týchto hodnôt. Výhodou (oproti programovaniu v C) knižnice Numba je, že po vyvolaní funkcie kernelu jej beh neskončí, pokým sa výpočet na GPU neukončí, a všetky údaje v pamäti nebudú synchronizované.

Python
import numpy
from numba import cuda


@cuda.jit
def my_kernel(io_array):
    """
    Code for kernel.
    """
    # code here
​

# Create the data array - usually initialized some other way
data = numpy.ones(256)

# Set the number of threads in a block
threadsperblock = 32 

# Calculate the number of thread blocks in the grid
blockspergrid = (data.size + (threadsperblock - 1)) // threadsperblock

# Now start the kernel
my_kernel[blockspergrid, threadsperblock](data)

# Print the result
print(data)

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

						import numpy
from numba import cuda
 
 
@cuda.jit
def my_kernel(io_array):
    """
    Code for kernel.
    """
    # code here
​
 
# Create the data array - usually initialized some other way
data = numpy.ones(256)
 
# Set the number of threads in a block
threadsperblock = 32 
 
# Calculate the number of thread blocks in the grid
blockspergrid = (data.size + (threadsperblock - 1)) // threadsperblock
 
# Now start the kernel
my_kernel[blockspergrid, threadsperblock](data)
 
# Print the result
print(data)

					

			

Veľkosť bloku vlákien

Dvojúrovňová hierarchia pracovných vlákien je dôležitá z dvoch dôvodov:

Na strane softvéru — veľkosť bloku určuje, koľko vlákien zdieľa oblasť zdieľanej pamäte.
Na strane hardvéru — veľkosť bloku určuje efektívne vyťaženie tzv. vykonávacej jednotky (nie všetky vlákna totiž môžu bežať súbežne; ich počet je obmedzený hardvérovými možnosťami zariadenia; podrobnosti viď v manuáli CUDA C Programming Guide).

Veľkosť bloku (t.j. počet vlákien tvoriacich blok) závisí na mnohých faktoroch, medzi ktoré patria:

Veľkosť poľa s údajmi,
veľkosť zdieľanej pamäte pre blok (toto je dané konkrétnym hardvérom, napr. 64kB),
maximálny počet vlákien podporovaných hardvérom v bloku (napr. 512 alebo 1024),
maximálny počet vlákien v rámci multiprocesora (MP), napr. 2048,
maximálny počet blokov v rámci MP, napr. 32,
počet vlákien, ktoré môžu bežať súbežne (tzv. „warp“), pre všetky doterajšie verzie CUDA je to hodnota 32.

Všetky vlákna vo warpe vykonávajú ten istý programový kód. Táto vlastnosť má obrovský dopad na efektivitu výpočtu. Ak totiž vykonávajú tú istú inštrukciu všetky vlákna, všetky bežia paralelne. Ale ak niektoré vlákno (dokonca viaceré) vykonávajú inú inštrukciu, warp sa musí rozdeliť na skupiny vlákien podľa vykonávaných inštrukcií, a inštrukcie týchto skupín sa vykonávajú sériovo, čo znižuje efektivitu výpočtu.

Niekoľko pravidiel pre odvodenie počtu vlákien v bloku:

Mal by to byť násobok hodnoty 32 (kvôli warpu),
na optimálne určenie hodnoty je vhodné vykonať niekoľko experimentov s rôznymi nastaveniami,
každý SP (streaming multiprocessor) na GPU by mal pre maximalizáciu svojho výkonu vyťažovať čo najviac warpov. Warp je totiž najmenšia plánovacia jednotka vlákien na GPU (pre systémy CUDA, ovšem). Veľkosť bloku by mala byť teda zvolená tak, aby sa pri danom kerneli vyťažila čo najväčšia časť hardvéru zariadenia. Pre viac detailov viď CUDA Occupancy Calculator spreadsheet.

Identifikácia vlákien

Telo kernel funkcie sa vykonáva práve toľkokrát, koľko vlákien bolo požadovaných konfiguráciou pri spustení kernelu. Preto je dôležité vedieť medzi sebou odlíšiť jednotlivé vlákna. Zvlášť, ak požadujeme, aby každé vlákno spracovalo inú časť údajov na vstupe.

Pre ľahšiu manipuláciu s viacrozmernými poliami umožňuje CUDA definovať viacrozmerné bloky a gridy. V príklade uvedenom vyššie je možné definovať premenné blockspergrid a threadsperblock ako n-tice pozostávajúce z jedného, dvoch alebo troch prirodzených čísiel. Z pohľadu efektivity vykonávania kódu je jedno, koľko rozmerné n-tice sa použijú na definíciu gridu a/alebo bloku. Avšak pri písaní programov, ktoré na vstupe narábajú s viacrozmernými poliami to môže byť dobrou pomôckou (jednak pre čitateľnosť takého kódu, na druhej strane aj pre zníženie rizika vnášania chýb pri prepočte indexov medzi rôznymi dimenziami).

Jednou z najjednoduchších možností identifikácie údajov, ktoré má dané vlákno spracovať (predpokladáme, že vstupné údaje sú v jednorozmernom poli), je výpočet indexu do poľa pomocou pozície vlákna v rámci bloku a bloku v rámci gridu:

Python
@cuda.jit
def my_kernel(io_array):
    # Thread id in a 1D block
    tx = cuda.threadIdx.x
    # Block id in a 1D grid
    ty = cuda.blockIdx.x
    # Block width, i.e. number of threads per block
    bw = cuda.blockDim.x
    # Compute index inside the array
    pos = tx + ty * bw
    # Check array boundaries
    if pos < io_array.size:
        # do the computation
        io_array[pos] *= 2

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14

						@cuda.jit
def my_kernel(io_array):
    # Thread id in a 1D block
    tx = cuda.threadIdx.x
    # Block id in a 1D grid
    ty = cuda.blockIdx.x
    # Block width, i.e. number of threads per block
    bw = cuda.blockDim.x
    # Compute index inside the array
    pos = tx + ty * bw
    # Check array boundaries
    if pos < io_array.size:
        # do the computation
        io_array[pos] *= 2

					

			

Poznámka: Ak veľkosť vstupného poľa údajov nie je násobkom bloku a gridu, vždy je nutné kontrolovať hranice, inak príde ku chybe prístupu k údajom mimo poľa!!!

Numba ponúka nasledovné „vstavané“ objekty pre zisťovanie geometrie hierarchie vlákien a pozíciu aktuálne vykonávaného vlákna v rámci tejto hierarchie:

numba.cuda.threadIdx: Index aktuálne vykonávaného vlákna v rámci bloku. Pre jednorozmerný blok je k dispozícii atribút x tohto objektu, pre dvoj rozmerný blok x a y a pre trojrozmerný x, y a z. Tieto atribúty nadobúdajú hodnoty od 0 po numba.cuda.blockDim-1.

numba.cuda.blockDim: Rozmery bloku dané pri volaní kernel funkcie. Táto hodnota je rovnaká pre všetky bloky gridu počas celej doby vykonávania kernel funkcie.
numba.cuda.blockIdx: Index bloku v gride. Podobne ako pri numba.cuda.threadIdx má tento objekt maximálne tri atribúty (x, y a z), v závislosti od počtu rozmerov gridu. Tieto atribúty nadobúdajú hodnoty od 0 po numba.cuda.gridDim-1.
numba.cuda.gridDim: Rozmery gridu dané pri volaní kernel funkcie.

Poznámka: Všetky tieto objekty môžu byť 1 až 3 rozmerné. Iba pre jednoduchosť zápisu sme vyššie miestami vynechali špecifikáciu rozmeru daného objektu.

Objekty blockDim, gridDim, threadIdx a blockIdx sú štandardnou súčasťou knižnice CUDA. Avšak Numba uľahčuje manipuláciu s geometriou a indexovaním pomocou ďalších nástrojov:

numba.cuda.grid(ndim): Vráti absolútnu pozíciu aktuálneho vlákna v celej mriežke blokov. Parameter ndim by sa mal zhodovať s počtom rozmerov definovaných pri volaní kernel funkcie. Ak je ndim 1, vráti sa iba jedno číslo. Ak má hodnotu 2 alebo 3, vráti sa príslušná n-tica čísel.

numba.cuda.gridsize(ndim): Vráti absolútnu veľkosť gridu (v počte vlákien). Znovu, aj tu platia rovnaké podmienky pre premennú „ndim“ ako pri predošlej funkcii.

Pomocou týchto funkcií je možné prepracovať (a výrazne zjednodušiť!) vyššie uvedený príklad nasledovne:

Python
@cuda.jit
def my_kernel2(io_array):
    pos = cuda.grid(1)
    if pos < io_array.size:
        # do the computation
        io_array[pos] *= 2

				
					
				1
2
3
4
5
6

						@cuda.jit
def my_kernel2(io_array):
    pos = cuda.grid(1)
    if pos < io_array.size:
        # do the computation
        io_array[pos] *= 2

					

			

Kompletný program vrátane vyvolania kernel funkcie by mal nasledovnú podobu:

Python
from __future__ import division
from numba import cuda
import numpy
import math


# CUDA kernel
@cuda.jit
def my_kernel(io_array):
    pos = cuda.grid(1)
    if pos < io_array.size:
        io_array[pos] *= 2 # do the computation


# Host code   
data = numpy.ones(256)
threadsperblock = 256
blockspergrid = math.ceil(data.shape[0] / threadsperblock)
my_kernel[blockspergrid, threadsperblock](data)
print(data)

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

						from __future__ import division
from numba import cuda
import numpy
import math
 
 
# CUDA kernel
@cuda.jit
def my_kernel(io_array):
    pos = cuda.grid(1)
    if pos < io_array.size:
        io_array[pos] *= 2 # do the computation
 
 
# Host code   
data = numpy.ones(256)
threadsperblock = 256
blockspergrid = math.ceil(data.shape[0] / threadsperblock)
my_kernel[blockspergrid, threadsperblock](data)
print(data)

					

			

Pred spustením programu nezabudnite nastaviť premennú NUMBA_ENABLE_CUDASIM na hodnotu 1, pokiaľ v systéme nemáte žiadnu grafickú kartu s podporou CUDA výpočtu!

Po dokončení behu programu by ste mali vidieť na výstupe:

[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.]

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

						[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  2.]

					

			

Modifikácia úlohy

Vieme, že pomocou x, y = cuda.grid(2) vieme získať indexy aktuálneho vlákna. Vašou úlohou je prepísať kernel funkciu tak, aby pracovala s dvoj-rozmerným poľom na vstupe. Pri prepisovaní kódu nezabudnite skontrolovať, či sú oba indexy v hraniciach poľa! Namiesto io_array.size použite io_array.shape. Vami napísanú kernel funkciu spustite pomocou nasledovného kódu:

Python
from __future__ import division
from numba import cuda
import numpy
import math


@cuda.jit
def my_kernel_2D(io_array):
    x, y = cuda.grid(2)
    ### YOUR SOLUTION HERE


data = numpy.ones((16, 16))
threadsperblock = (16, 16)
blockspergrid_x = math.ceil(data.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(data.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)
my_kernel_2D[blockspergrid, threadsperblock](data)
print(data)

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

						from __future__ import division
from numba import cuda
import numpy
import math
 
 
@cuda.jit
def my_kernel_2D(io_array):
    x, y = cuda.grid(2)
    ### YOUR SOLUTION HERE
 
 
data = numpy.ones((16, 16))
threadsperblock = (16, 16)
blockspergrid_x = math.ceil(data.shape[0] / threadsperblock[0])
blockspergrid_y = math.ceil(data.shape[1] / threadsperblock[1])
blockspergrid = (blockspergrid_x, blockspergrid_y)
my_kernel_2D[blockspergrid, threadsperblock](data)
print(data)

					

			

Na výstupe by ste mali uvidieť:

Python
[[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]]

				
					
				1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

						[[ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]
 [ 2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.  2.]]

					

			

Násobenie matíc

Záujemcom o zložitejší príklad CUDA výpočtu odporúčame pokračovať sekciou „A more complex example: matrix multiplication“ na stránke https://nyu-cds.github.io/python-numba/05-cuda/

Pre pochopenie alokácie pamäte v programovom kóde hosta je vhodné mrknúť na uvedenej stránke aj na odsek vyššie: „Memory management“.