W XX wieku procesory graficzne służyły do transformacji grafiki, głównie dwuwymiarowej. Potrzeby przemysłu związanego z nowoczesnymi grami wykorzystującymi grafikę trójwymiarową doprowadziły do intensywnego rozwoju i co równie ważne - upowszechnienia się, wyspecjalizowanych procesorów zwanych GPU (Graphics Processor Units), zdolnych renderować coraz bardziej skomplikowane sceny. W 2001 wraz z powstaniem technologii Pixel Shader i Vertex Shader zaczęto wykorzystywać procesor graficzny do obliczeń niezwiązanych z grafiką, określane mianem GPGPU - General-Purpose computing on Graphics Processor Units.
W 2004, firma Ageia zaprezentowała przełomowy produkt, znany pod nazwą PhysX a będący sprzętowym akceleratorem obliczeń związanych z symulacjami fizyki w grach komputerowych. Po co fizyka w grach? Rozwój grafiki 3D szybko doprowadził do sytuacji w której komputer mógł całkiem rozsądnie wyrenderować żądaną scenę. Jednak aby to zrobił musi wiedzieć co renderować. Po części informacje o położeniu objektów w grach są zapisywane z odgrywanych przez aktorów scen lub wprowadzane ręcznie przez animatora. Jednak ma to wiele wad: jest kosztowne, zabiera duzo pamięci a powtarzalność ruchów negatywnie wpływa na odbiór akcji. Dlatego optymalnym rozwiązaniem było generowanie scen przez pewne algorytmy. I tu przeszkodą stał się nasz mózg, który jest siecią neuronową doskonale znającą podstawowe prawa fizyki! Bez trudu potrafimy odgadnąć jak ruszał by się słoń, gdyby był ze styropianu, znamy prawa odbicia, wiemy jak przepływa ciecz itp. Nie łatwo nas oszukać - więć programiści gier komputerowych postanowili nie "generować", ale (prawie uczciwie) symulować świat w grach. Jednak, aby obliczyć trajektorie na przykład miliona liści opadających z drzewa, potrzeba ogromnych mocy obliczeniowych. Z drugiej strony gra wymaga wykonania tych obliczeń w pewnym ogranicznonym czasie. Właśnie ta potrzeba, wraz z siedemdziesięcio-miliardowym biznesem gier, doprowadziła opracowania koncepcji sprzętowego akceleratora fizyki- tzw. PPU (Physics Processing Unit).
Co się stało, że firmy Ageia nie ma już na runku? Okazało się, że producenci kart graficznych z procesorami GPU, zauważyli, że jest to ich branża. Ponadto, okazało się, że architektura GPU jest we wielu aspektach zliżona do architektury PPU. GPU wykonuje równoległe obliczenia na "pixelach", PPU wykonuje pewne obliczenia na wielu cząstkach. Firma Nvidia wykonała w 2006 kilka kroków - po pierwsze udostępniła standart programowania jednostek GPU zwany CUDA i zachęciła świat nauki to rozwijania oprogramowania. Po drugie wykupiła Agei-ę i zaimplementowała bibliotekę PhysX w CUDA, wobec czego jej akceleracja była możliwa nie tylko z procesorami PPU, ale na każdej nowej karcie firmy Nvidia. Praktycznie od momentu pojawienia się generacji GeForce serii 8, wszystkie procesory GPU z Nvidii są kompatybilne z CUDA i mogą wykonywać równolegle napisane programy. Różnice pomiędzy procesorami polegają głównie na ilości rdzeni, zwanych "CUDA -cores", których liczba waha się od 16 do 2500!
Moc obliczeniowa procesorów GPU a także, co jest niesłychanie ważne - przepustowość pamięci są o wiele większe od parametrów dostępnych na CPU. Oczywiście jest to kosztem ograniczenia możliwości procesorów GPU. Aby zdać sobie sprawę z liczb, warto przyjrzeć się grafikom opublikowanym przez Nvidię:
Jak już się dowiedzieliśmy, nowoczesne procesory graficzne to urządzenia obliczeniowe, zdolne do wykonywania trylionów operacji zmiennoprzecinkowych na sekundę.
Na dzień dzisiejszy dostępne urządzenia GPU działające w standardzie CUDA można podzielić na trzy generacje o nazwach kodowych: Tesla, Fermi i Kepler. Każda z nich oferuje coraz bardziej zaawansowane funkcje i lepszą wydajność.
Procesor GPU jest zorganizowany wokół koncepcji multiprocesorów strumieniowych (MP). Takie multiprocesorowy składają się z kilku-kilkunastu procesorów skalarnych (SP), z których każdy jest zdolny do wykonywania wątku w SIMT (Single Instruction, Multiple Threads). Każdy MP posiada również ograniczoną ilość wyspecjalizowanej pamięci on-chip: zestaw rejestrów 32-bitowych, wspólny blok pamięci i pamięci podręcznej L1, cache stałych i cache tekstur. Rejestry są logicznie lokalne dla procesora skalarnego, ale inne typy pamięci są wspólne dla wszystkich SPs w każdym w MP. Umożliwia to m.in. wymianę danych pomiędzy wątkami operującymi na tym samym MP.
Procesory graficzne mogą różnić się ilością i organizacją rdzeni obliczeniowych (CUDA-cores). Implementacja algorytmów w równoległy sposób powinienna być niezależna od konkrentego sprzętu na którym się one wykonują. W ogólności jest to zagadnienie - jak się wydaje - nierozwiązywalne, ale technologia CUDA uczyniła olbrzymi krok w kierunku uzyskania niezależności od sprzętu.
Po pierwsze programując na GPU tworzymy dużo więcej wątków niż jest fizycznie dostępnych rdzeni obliczeniowych. System kolejkowania wykonań zajmuje się częściową serializacją wywołań, tym intensywniejszą im gorszy sprzęt posiadamy.
Po drugie dostęp do sprzętu jest realizowany przez warstwe abstrakcji, zwaną właśnie CUDA, która zawiera w sobie informacje o hierarchii pamięci. Programista, nie przejmuje się tym ile rdzeni wykonuje jego program, ale wie, że pewne grupy wątków mają dostęp do wspólnej szybkiej pamięci.
Programy napisane na CUDA dramatycznie różnią się od wielowątkowych równoległych odpowiedników napisanych na nawet kilkunastordzeniowe jednostki CPU. Z jednej strony jest to ograniczenie, skutkujące koniecznością przepisania wszystkich algorytmów. Jednak poztywnym efektem programowania na takiej architekturze jest ogromna kompatybilność i przenośność kodu na współczesne i z dużym prawdopodobieństwem przyszłe urządzenia.
Być może najbardziej istotną cechą architektury CUDA jest hierarchia pamięci z różnicą czasów dostępu o 1-2 rzędy wielkości pomiędzy kolejnymi poziomami. Najwolniejszą z punktu widzenia GPU jest pamięć RAM komputera gospodarza (host). Pamięć ta jest oddzielona od procesora graficznego magistralą PCIe z teoretyczną maksymalną przepustowością w jednym kierunku 16 GB/s (PCI Express 3.0,x16).
Następną w kolejności jest globalna pamięć urządzenia GPU, który jest obecnie ograniczony do kilku gigabajtów (najnowsze karty mają juz 12GB) i o szerokości pasma około 100-200 ~ Gb/s. Jest to bardzo duża wartość, jednak dostęp do globalnej pamięci jest operacją wysokiej latencji wynoszącą kilkaset cykli zegara GPU.
Najszybszym rodzajem pamięci dostępnym obecnie na GPU jest pamięć współdzielona (shared memory) znajdująca się bezpośrednio na multiprocesorze. Jest ona obecnie ograniczona do 48 kB, ale ma przepustowość około ~1,3 TB/s. Co ciekawsze, latenacja w dostępie do tej pamięci jest bardzo mała - podobna jak dostęp do rejestrów jednostek SP.
Powyższy opis łatwo sugeruje strategię pisania wydajnych programów na CUDA, a które można streścić jako: przenieść jak najwięcej danych, jak to możliwe do najszybszych rodzaju dostępnej pamięci i przechowywać je tam tak długo, jak to możliwe, przy jednoczesnej minimalizacji dostępu do wolniejszych rodzajów pamięci. Dodatkowo, jeśli dostęp do wolniejszej pamięci jest konieczny, można próbować wykorzystać "wolny" czas na wykonanie operacji arytmetycznych na pozostałych wątkach.
W programowaniu na CUDA czy OpenCL wykorzystujemy wątki. Co jest niespotykane w pisaniu kodu HPC na zwykłe procesory, wątków powinno się włączać o wiele więcej od dostępnych fizycznych rdzeni procesora. W praktyce programowania na CPU znany jest fakt, że nadmierne rozmnożenie wątków zazwyczaj spowalnia wykonanie programu, ze względu na "context switching", który może prowadzić do nieefektywnego wykorzystania pamieci cache jednostki centralnej a ponadto prowadzi w nieunikniony sposób do zmniejszenia dostępnej pamięci RAM. W programowaniu w standardzie CUDA wątki są lekkie, ich przełączanie nie zabiera więcej niż pojedyńczych cykli procesora. Co więcej, nadmiarowe wątki potrafią mieć pozytywny efekt na wydajność, gdyż mogą przesłonić latencję dostępu do pamięci głównej.
Z punktu widzenia użytkownika, programy CUDA są zorganizowane w jądrach. Jądro Jest to funkcja, która jest wykonywana wielokrotnie, jednocześnie na różnych multiprocesorach. Każda instancja jądra zwana jest wątkiem i jest przydzielana do pojedyńczego procesora skalarnego (SP). Wątki są grupowane w jedno-, dwu - lub trójwymiarowe bloki przypisanych do multiprocesorów jeden do jednego - czyli mamy gwarancję, że w ramach jednego bloku jesteśmy na tym samym MP.
Kernele pisze się w okrojonym dialekcie C/C++, zwanym CUDA-C. Referencyjnym i stosunkowo przystępnie napisanym dokumentem jest CUDA C Programming Guide, który jest dostępny pod adresem: http://docs.nvidia.com/cuda/cuda-c-programming-guide/ W tym kursie nie będziemy systematycznie analizować wszystkich elementów CUDA-C. Postąpimy inaczej. Pokażemy na przykładach typowe zastosowania programowania GPU, które są niezwykle przydatne w fizyce. Rozwiązanie własnego problemu, proponujemy rozpocząć od znalezienia stosownego przykładu i próbie samodzielnego dostosowania jego kodu. Jest to rozwiązanie typu "crash course", które nie zastąpi systematycznego kursu. Jednak w dużej większości przypadków, doświadczenie nabyte podczas zabawy z zamieszczonymi przykładami umożliwi optymalne rozwiązanie napotkanego problemu.
W tych materiałach zostanie zaprezentowane podejście, które pozwoli maksymalnie uprościć drogę do efektywnej pracy w CUDA. W tym celu zostanie zastosowany pakiet pyCUDA, który umożliwi pracę interaktywną z urządzeniem CUDA bez kompromisu jeśli chodzi o wydajność. W pyCUDA, sposób pracy sprowadza się do napisania jądra obliczeniowego w języku C z rozszerzeniami CUDA. Wywołanie jądra, jego kompilacja oraz inspekcja danych, w tym transfer do i z urzadzenia, wykonuje się w wygodnym języku python.
Przedstawimy teraz pierwszy program, który będzie napisany w CUDA, z użyciem pyCUDA. Przykład ten będzie bardzo prosty: mnożenie wektora przez liczbę. Na tym prostym przykładzie poznamy sposób pracy, który potem będziemy mogli użyć jako szablonu do tworzenia bardziej zaawansowanych programów.
Najprostszy sposób inicjalizacji urządzenia GPU tak by móc go dalej używać wygląda następująco:
In [1]:
import numpy as np
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
Oprócz modułow z pakietu pyCUDA, importujemy również numpy. Jest to niezwykle istotne, bowiem tablice numpy to podstawowa struktura danych którą będziemy przesyłać na urządzenie GPU i z powrotem.
Utwórzmy wektor danych, na przykład wykorzystując linspace:
In [2]:
a = np.linspace(1,16,16).astype(np.float32)
a
Out[2]:
Celem naszego programu będzie pomnożenie wszystkich elementów tej tablicy przez zadaną liczbę.
Tablica a znajduje się w pamięci systemu gospodarza. Aby wykonać operacje na niej z pomocą GPU musimy zaalokować pamięć na urządzeniu GPU:
In [3]:
a_gpu = cuda.mem_alloc(a.nbytes)
In [ ]:
a następnie skopiować dane do GPU:
In [4]:
cuda.memcpy_htod(a_gpu, a)
Mamy teraz na GPU zainicjalizowaną tablice, do ktorej możemy się z poziomu hosta (gospodarza) odwoływać przez a_gpu. Typową praktyką w programowaniu na GPU jest posiadania dwóch kopii danych - jednej w pamięci dostępnej dla CPU a drugą w pamięci urządzenia GPU. Warto jeszcze zaznaczyć, że nasza tablica a_gpu znajduje się w pamięci globlanej GPU.
Teraz musimy napisać jądro, które będziemy wywoływać na urządzeniu na danych w pamięci globalnej GPU. Jądro jest napisane w zmodyfikowanej wersji języka C, zwanej CUDA-C.
Jądro, które będzie mnożyło wszystkie elementy tablicy przez 2 ma następującą postać:
__global__ void dubluj(float *a)
{
int idx = threadIdx.x;
a[idx] *= 2;
}
__global__ jest rozszerzeniem CUDA C.a_gpu a każda będzie mnożyła jeden element wektora.threadIdx jest strukturą określającą numer wątku wewnątrz bloku, który wykonuje daną kopię. Oprócz niego ważny jest też blockIdx, ale w tym przypadku odpalamy tylko jeden blok więc nie jest nam potrzebny. Indeksy wątku, czy też wewnątrz bloku czy numer bloku, to jedyny sposób rozróżnienia wątków! Zmienna idx jest liniowym wskaźnikiem, który będzie się zmieniał od 0 do N-1.W pyCUDA źródło modułu podajemy jako string funkcji SourceModule w poniższy sposób:
In [38]:
mod = SourceModule("""
__global__ void dubluj(float *a)
{
int idx = threadIdx.y;
a[idx] *= 3.12;
}
""")
Zmienna mod jest teraz obiektem, który ma m.in. funkcję mod.get_function("nazwa") zwracającą funkcję pythonową, której wywołanie uruchomi odpowiednie jądro ze źródła uprzednio podanego do SourceModule. Sprawdźmy:
In [39]:
func = mod.get_function("dubluj")
Wywołajmy teraz jądro!
In [34]:
func(a_gpu, block=(16,1,1))
W tablicy a_gpu powinniśmy mieć teraz
In [35]:
print( a)
cuda.memcpy_dtoh(a, a_gpu)
print( a)
Powyższy sposób budowania programu w CUDA jest bardzo podobny gdybyśmy używali jedynie kompilatora oraz języka C/C++ zamiast pythona z pyCUDA.
Warto zaznajomić się z kilkoma udogodnieniami, które mamy wbudowane w pyCUDA. Pierwszym z nich jest zautomatyzowanie procesu alokacji i kopiowania danych. Jeśli chcemy wykonać po kolei: alokacje wektora na GPU, transfer danych z wektora numpy, wywołanie jądra oraz nadpisanie wyjściowego wektora numpy wynikiem z GPU, to możemy użyć funkcji InOut:
In [41]:
func(cuda.InOut(a), block=(4, 4, 1))
print( a)
gpuarrayModuł gpuarray umożliwia bardzo zwarty zapis operacji na wektorach na GPU. Kluczową komendą jest gpuarray.to_gpu(.., która pozwala skopiować wektor numpy do urządzenia GPU, zwracając uchwyt do danych na GPU. Pewne operacje można wykonać automatycznie na GPU po prostu wpisując wzór arytmetyczny.
In [11]:
import pycuda.gpuarray as gpuarray
a = np.linspace(1,16,16).astype(np.float32)
a_gpu = gpuarray.to_gpu(a.astype(np.float32))
a_doubled = (a_gpu*2).get()
print (a_gpu)
print (a_doubled)
Uwaga printf - działa tylko z konsoli!
In [12]:
%%writefile cuda_printf.py
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
mod = SourceModule("""
#include <stdio.h>
__global__ void printf_test()
{
printf("GONZO %d.%d.%d\\n", threadIdx.x, threadIdx.y, threadIdx.z);
}
""")
func = mod.get_function("printf_test")
func(block=(2,2,1))
In [23]:
W CUDA wątki można uruchamiać w blokach (blocks). Należy zapamiętać regułe:
wszystkie wątki na jednym bloku zawsze będą uruchamiane na tym samym multiprocesorze
Co za tym idzie, będą miały do dyspozycji tą samą pamięc typu "shared memory". Maksymalna liczba wątków w jednym bloku jest z reguły ograniczona do 512 lub 1024 w zależności od typu urządzenia GPU. Warto zauważyć, że jest to liczba zdecydowanie większa od liczby procesorów skalarnych w na jednym multiprocesorze.
Jeżeli chcemy odpalic np. milion wątków, to musimy wykorzystać tzw. grid bloków. Czyli odpalając jądro obliczeniowe specyfikujemy ile potrzeujemy wątków i ile jaki chcemy mieć rozmiar bloku. Chcąc mieć np. $128\times 64$ wątki, rozsądnym jest odpalenie jądra z rozmiarem gridu $128$ i rozmiarem bloku $64$.
Ponadto, ponieważ często operujemy siatkami reprezentujemy dwu i trój-wymiarowe pola, CUDA ma wbudowany mechanizm, który ułatwia posługiwaniem się dwu i trój-wymiarowymi siatkami.
Poniższy przykład pokazuje jak zbudować grid składający się z dwóch bloków, po pięc wątków.
Zastąp w kodzie if(threadIdx.x==1) przez:
if(blockIdx.x==1)if(idx==1)i sprawdź działanie.
In [45]:
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
mod = SourceModule("""
__global__ void kernel(float *a)
{
int idx = threadIdx.x + blockDim.x*blockIdx.x;
if(threadIdx.x==0)
a[idx] += 1.0f;
a[idx] += 10.0f * idx;
}
""")
a = np.zeros(10).astype(np.float32)
func = mod.get_function("kernel")
print (np.linspace(0,9,10))
print ("----------------")
print (a)
func(cuda.InOut(a),block=(5,1,1),grid=(2,1,1))
print (a)