Programowanie równoległe w języku CUDA C - #5 - Biblioteka CUBLAS

[Qitr 41]

CUDA na kiju (a nawet GPU) cz 5.

Wstęp

W poprzednich częściach naszego cyklu CUDA na kiju działaliśmy bardzo niskopoziomowo tzn. sami tworzyliśmy funkcję jądra definiowaliśmy wymiary bloku itd. Czy zawsze tak trzeba robić? Nie, wszyscy, którzy mieli jakiekolwiek wykłady związane z programowaniem wiedzą (a ci którzy nie mieli zaraz się dowiedzą), że dobry programista to leniwy programista. Zatem jeśli istnieją gotowe biblioteki zawierające funkcję, która pozwala na rozwiązanie postawionego nam problemu to z reguły nie napiszemy tego kodu lepiej niż grupa zawodowych informatyków. Dlatego też w tym artykule sprawdzimy jak z mnożeniem macierzy o wymiarach 1021x1021 poradzi sobie biblioteka CUBLAS.

Ten artykuł bierze udział w naszym konkursie! 
Na zwycięzców czekają karty podarunkowe Allegro, m.in.: 2000 zł, 1000 zł i 500 zł.

Potrafisz napisać podobny poradnik? Opublikuj go na forum i zgłoś się do konkursu!
Czekamy na ciekawe teksty związane z elektroniką i programowaniem. Sprawdź szczegóły »

CUBLAS

Nazwa biblioteki pochodzi od CUDA BLAS (Basic Linear Algebra Subropgrams), co oznacza, że jest to biblioteka do algebry liniowej działająca z wykorzystaniem języka CUDA C, a zatem również procesora graficznego. Ze względu na zgodność CUBLAS z jego pierwowzorem (BLAS) przechowuje on macierze kolumnami i indeksuje je od 1 (!!!), a nie od 0 jak ma to miejsce w języku C. Niby drobny szczegół, a jednak trzeba uważać. Dodatkowo należy pamiętać, że procedury działają „w miejscu” tzn. nadpisują swoje wyniki tam gdzie to tylko możliwe np. kombinacja liniowa wektorów x i y nadpisze wartość wektora y:

 (1)

Aczkolwiek w przypadku mnożenia macierzy nie zawsze jest to możliwe, gdyż wymiary macierzy wynikowej mogą być inne niż mnożonych macierzy. Nie zmienia to faktu, że często jest to dobry pomysł aby nadpisywać wyniki, gdyż minimalizuje to liczbę transferów danych od gospodarza do urządzenia i vice versa.

W BLAS występują trzy poziomy abstrakcji algorytmów:

Na pierwszym poziomie znajdują się operacje na dwóch wektorach. Ilość danych oraz obliczeń wzrasta liniowo wraz z rozmiarem problemu, ale CMGA jest mniej więcej stałe. Na drugim poziomie znajdują się operacje na macierzy i wektorze. Ilość danych i obliczeń wzrasta kwadratowo wraz z rozmiarem problemu, ale CGMA również nie zależy od rozmiaru problemu. Na trzecim poziomie znajdują się operacje na macierzach. Ilość danych wzrasta kwadratowo wraz z rozmiarem problemu, ilość obliczeń wzrasta proporcjonalnie do trzeciej potęgi rozmiaru problemu, a CGMA wzrasta liniowo wraz z rozmiarem problemu. Z części drugiej naszego cyklu wiemy, że aby opłacało się zatrudniać GPU do pracy CGMA musi być odpowiednio dużo. Z tego względu największy zysk z zastosowania CUBLAS-a będzie jeśli użyjemy go do rozwiązania problemu z trzeciego poziomu. Zysk ten zwiększymy jeśli zastosujemy procedury biblioteczne dla całego bloku macierzy (np. pomnożymy 1000 macierzy na raz).

Konwencja nazw

Nazwy procedur BLAS pochodzą z zamierzchłych czasów kiedy to identyfikatory FORTRAN-u miały maksymalnie 6 znaków, z których pierwszy oznaczał typ danych. W bibliotece CUBLAS właściwą nazwę procedury (pisaną z wielkiej litery – camelCase) poprzedza się przedrostkiem cublas. Na początku nazwy występuje prefiks t oznaczający typ danych:

• H - __half – typ połówkowej precyzji, gdy mamy mało pamięci
• S – float – typ pojedynczej precyzji (silna strona CUDA)
• D – double – typ podwójnej precyzji
• C – cuComplex – typ zespolony pojedynczej precyzji
• Z – cuDoubleComplex – typ zespolony podwójnej precyzji

Dla procedur poziomu drugiego i trzeciego (a nimi będziemy się tutaj zajmować) po prefiksie typu danych t mamy rodzaj macierzy występującej w obliczeniach:

• tGE – General – dowolna macierz
• tGB – General Band – dowolna macierz pasmowa
• tSY – Symmetric – macierz symetryczna
• tSP – Symmetric Packed – symetryczna spakowana
• tSB – Symmetric Band – symetryczna pasmowa
• tHE – Hermitian – macierz Hermitowska
• tHP – Hermitian Packed - macierz Hermitowska spakowana
• tHB – Hermitian Band - macierz Hermitowska pasmowa
• tTR – Triangular – macierz trójkątna
• tTP – Triangular Packed – macierz trójkątna spakowana
• tTB – Triangular Band – macierz trójkątna pasmowa
• tDG – DiaGonal – macierz diagonalna

Jak widać trzeba łyknąć nieco matematyki, żeby w ogóle skumać co mamy do wyboru. W naszym przypadku zajmować się będziemy mnożeniem dwóch dowolnych (bo w końcu wygenerowanych losowo) macierzy wypełnionych danymi pojedynczej precyzji. Stąd przedrostek cublasSge.

Procedury poziomu 3

W opisie procedur wykorzystam operator ◊, który może mieć następujące znaczenie:

• Tożsamości  – CUBLAS_OP_N
• Transpozycji  – CUBLAS_OP_T
• Sprzężonej transpozycji  – CUBLAS_OP_C

Spis wszystkich procedur poziomu trzeciego umieszczam poniżej:

Nas interesuje mnożenie macierzy czyli Matrix-Matrix product ze stałymi α=1 i β=0. Tak więc podsumowując ten dosyć rozwlekły aczkolwiek wydaje się niezbędny wstęp teoretyczny wykorzystamy funkcję: cublasSgemm(), która przyjmuje następujące argumenty:

• cublasHandle_t handle – uchwyt do biblioteki CUBLAS.
• cublasOperation_t transa – czyli informacja o tym czy macierz A ma być transponowana/transponowana ze sprzężeniem – umieszczamy tutaj informacje o wartości operatora ◊ dla macierzy A
• cublasOperation_t transb – czyli informacja o tym czy macierz A ma być transponowana/transponowana ze sprzężeniem – umieszczamy tutaj informacje o wartości operatora ◊ dla macierzy B
• int m – liczba wierszy macierzy A po zastosowaniu operatora ◊, która jest równa liczbie wierszy macierzy C (co wynika z własności mnożenia macierzy)
• int n – liczba kolumn macierzy B po zastosowaniu operatora ◊, która jest równa liczbie kolumn macierzy C (co wynika z własności mnożenia macierzy)
• int k – liczba kolumn macierzy A po zastosowaniu operatora ◊, która jest równa liczbie wierszy macierzy B po zastosowaniu operatora ◊ (nie musze już chyba dodawać, że wynika to z własności mnożenia macierzy )
• const float * alpha -  skalar przez który mnożona jest wynik iloczynu A*B. Zwracam uwagę, że przekazujemy go przez wskaźnik! Czyli musimy zdefiniować zmienna alpha = 1.0 i przekazać do niej wskaźnik a nie wpisać po prostu 1.0 jako argument funkcji.
• const float * A – macierz A w pamięci urządzenia
• int lda – główny kierunek (głupio to po polsku brzmi – od angielskiego leading direction) dwu wymiarowej tablicy służącej do przechowywania macierzy A w pamięci. Podajemy liczbę kolumn w macierzy A
• const float * B – macierz B w pamięci urządzenia
• int ldb – główny kierunek dwu wymiarowej tablicy służącej do przechowywania macierzy B w pamięci. Podajemy liczbę kolumn w macierzy B
• const float * beta -  skalar przez który mnożona jest macierz C. W naszym przypadku ma on wartość zero. Zwracam uwagę, że przekazujemy go przez wskaźnik! Czyli musimy zdefiniować zmienna beta = 0.0 i przekazać do niej wskaźnik a nie wpisać po prostu 0 jako argument funkcji.
• const float * C – macierz C w pamięci urządzenia
• int ldc – główny kierunek dwu wymiarowej tablicy służącej do przechowywania macierzy C w pamięci. Podajemy liczbę kolumn w macierzy C

Kod programu

Po dość obszernym opisie biblioteki CUBLAS czas łyknąć nieco praktyki. Zajrzyjmy zatem do kodu:

1.	// Mnożenie macierzy
2.	// wersja GPU wykorzystująca bibliotekę CUBLAS
3.	 
4.	#include "cuda_runtime.h"
5.	#include "device_launch_parameters.h"
6.	#include "helper_cuda.h"//-I$(NVCUDASAMPLES_ROOT)/common/inc
7.	 
8.	#include "matmul.h"
9.	 
10.	#include <cublas_v2.h> // wymaga dołączenia biblioteki cublas.lib linker->input->dopisac cublas.lib
11.	// definiujemy makro ktore sprawdza status bledu operacji wykonanych przez cublas
12.	#define CUBLAS_SAFE_CALL(code)                                      \
13.	{                                                                   \
14.	  cublasStatus_t err = code;                                        \
15.	  if(CUBLAS_STATUS_SUCCESS != err)                                  \
16.	  {                                                                 \
17.	    fprintf(stderr, "CUBLAS error in file '%s' in line %i : %d.\n", \
18.	            __FILE__, __LINE__, err);                               \
19.	    exit(EXIT_FAILURE);                                             \
20.	  }                                                                 \
21.	}
22.	 
23.	void matmul(float *C, float *A, float *B,	int m, int p, int n)
24.	{
25.		checkCudaErrors(cudaSetDevice(0));//wybor urzadzenia
26.	 
27.		float *dev_A, *dev_B, *dev_C; //deklaracja macierzy w pamieci urzadzenia
28.		//Alokacja pamieci
29.		checkCudaErrors(cudaMalloc(&dev_A, m*p * sizeof(float)));
30.		checkCudaErrors(cudaMalloc(&dev_B, p*n * sizeof(float)));
31.		checkCudaErrors(cudaMalloc(&dev_C, m*n * sizeof(float)));
32.		//kopiowanie danych z pamieci gospodarza do pamieci urzadzenia
33.		checkCudaErrors(cudaMemcpy(dev_A, A, m*p * sizeof(float),cudaMemcpyHostToDevice));
34.		checkCudaErrors(cudaMemcpy(dev_B, B, p*n * sizeof(float),cudaMemcpyHostToDevice));
35.		//Nowosc - cublas
36.		cublasHandle_t h;//deklracja uchwytu do biblioteki
37.		CUBLAS_SAFE_CALL(cublasCreate(&h));//tworzymy uchwyt
38.		// definicja skalarow do mnozenia macierzy C<-alpha*AB + beta*C
39.		const float alpha = 1.0f;
40.		const float beta = 0.0f;
41.	 
42.		cudaEvent_t start, stop; // pomiar czasu wykonania jądra
43.		checkCudaErrors(cudaEventCreate(&start));
44.		checkCudaErrors(cudaEventCreate(&stop));
45.		checkCudaErrors(cudaEventRecord(start, 0));
46.	 
47.		/* Biblioteka CUBLAS przechowuje macierze kolumnami,
48.		   a nie wierszami, więc zmieniamy kolejność czynników
49.		   iloczynu macierzowego. */
50.		CUBLAS_SAFE_CALL(cublasSgemm(h, CUBLAS_OP_N, CUBLAS_OP_N,n, m, p, &alpha, dev_B, n, dev_A, p, &beta, dev_C, n));
51.	 
52.		checkCudaErrors(cudaGetLastError());// sprawdzmy status bledu
53.	 
54.		checkCudaErrors(cudaEventRecord(stop, 0));
55.		checkCudaErrors(cudaEventSynchronize(stop));
56.		float elapsedTime;
57.		checkCudaErrors(cudaEventElapsedTime(&elapsedTime,start, stop)); //obliczamy czas trwania jadra
58.		checkCudaErrors(cudaEventDestroy(start));
59.		checkCudaErrors(cudaEventDestroy(stop));
60.	 
61.		checkCudaErrors(cudaDeviceSynchronize());
62.		//Kopiowanie danych z pamieci urzadzenia do pamieci gospodarza
63.		checkCudaErrors(cudaMemcpy(C, dev_C, m*n * sizeof(float),cudaMemcpyDeviceToHost));
64.	 
65.		CUBLAS_SAFE_CALL(cublasDestroy(h));//niszczymy uchwyt do biblioteki
66.		//Zwalniamy zaalokowana na urzadzeniu pamiec
67.		checkCudaErrors(cudaFree(dev_C));
68.		checkCudaErrors(cudaFree(dev_B));
69.		checkCudaErrors(cudaFree(dev_A));
70.	 
71.		checkCudaErrors(cudaDeviceReset()); // dla debuggera
72.		//Wypisanie czasu trwania jadra
73.		printf("Czas wykonania jadra = %.3f ms (%6.3f GFLOP/s)\n",
74.			elapsedTime, 2e-6 * m * n * p / elapsedTime);
75.	}

Zanim skompilujemy program musimy pamiętać o dołączeniu biblioteki cublas.lib. W tym celu we właściwościach projektu wybieramy linker->input->additional dependences i tam dopisujemy cublas.lib:

Rysunek 1. Dołączanie biblioteki CUBLAS do projektu w VS 2017.

Teraz możemy spokojnie wykorzystać dyrektywę #include <cublas_v2.h>. Następnie definiujemy sobie makro sprawdzające status błędu w operacjach wykonywanych przez bibliotekę. W funkcji matmul() standardowo zaczynamy od wyboru urządzenia, alokacji pamięci i kopiowania danych do pamięci urządzenia. Potem pojawia się nowość. Deklaracja uchwytu do biblioteki oraz jego utworzenie (linie 36 i 37). Zawsze musimy wykonać te operacje kiedy pracujemy z biblioteką CUBLAS. Następnie definiujemy sobie skalary, potrzebne do funkcji cublasSgemm(). Potem standardowy mechanizm do pomiaru czasu trwania jądra, ale zamiast jawnego wywołania funkcji jądra wywołujemy wyżej już wspomnianą funkcję cublasSgemm() i tutaj pułapka! Musimy zamienić kolejność mnożenia w iloczynie macierzowym ze względu na fakt, że CUBLAS przechowuje macierze kolumnami (tak jak MATLAB), a język C wierszami. Dlaczego? Już spieszę z wyjaśnieniem: w MATLABIE wykonaliśmy operacje:

 (2)

A potem transponowaliśmy macierze ze względu na na to, że w języku C macierze przechowujemy wierszami, a w MATLAB-ie kolumnami. Jako, że CUBLAS przechowuje macierze tak samo jak MATLAB to aby otrzymać poprawny (z punktu widzenia sprawdzania w naszym CUDownym programie) wynik musimy obliczyć macierz C^T, którą można opisać w następujący sposób:

 (3)

Zatem skoro mamy już przetransponowane macierze A i B w naszym programie to prosta zamiana kolejności czynników załatwi sprawę .

Gdy już wykonamy mnożenie macierzy kopiujemy dane z pamięci urządzenia do pamięci gospodarza, a następnie niszczymy uchwyt do biblioteki CUBLAS. Potem standardowo zwalniamy pamięć urządzenia i wypisujemy na konsolę czas wykonania jądra.

Wyniki

Po zapoznaniu się z kodem nie pozostaje nam nic innego jak skompilować i sprawdzić efekty. Wciskamy zatem F5 i obserwujemy co się wydarzy. Jeśli wszystko wykonaliśmy poprawnie to naszym oczom powinna ukazać się konsola z następującym napisem:

Rysunek 2. Wyniki i czas działania biblioteki CUBLAS.

Porównajmy to sobie z wynikami naszej „poprawionej” wersji programu:

Rysunek 3. Wyniki i czas działania naszego kodu.

Jak widać jeśli chodzi o czas trwania jądra biblioteka CUBLAS jest bezkonkurencyjna. Totalnie nas zmiażdżyła. Czas wykonania mnożenia macierzy przy jej użyciu jest około 6 razy krótszy niż z wykorzystaniem naszej funkcji jądra. Jest jednak jedno ale – wywołanie biblioteki ma swój koszt. Czas całkowity trwania programu jest 1.5 raza dłuższy. Jaka z tego płynie dla nas nauka? Jeśli chcemy pomnożyć jedną macierz lepiej jest skorzystać z naszego kodu, który pomimo, że dłużej oblicza macierz wynikową to jest generuje zdecydowanie mniejszy narzut czasowy związany z wywaleniem funkcji jądra, ale jeśli chcemy pomnożyć np. 1000 macierzy wówczas zysk z zastosowania biblioteki CUBLAS będzie istotny.

Zakończenie

To już ostatni artykuł z naszego cyklu CUDA na kiju (a nawet GPU). Mam nadzieję, że wśród czytelników znajdą się osoby, którym moje wypociny przypadły do gustu . Opisane zagadnienia są dość podstawowe nie mniej jednak powinny dać dobre podstawy do dalszego rozwoju w tej dziedzinie. Dziękuję za uwagę. Piotrek się odmeldowuje

Spis treści:

W cyklu artykułów CUDA na kiju (a nawet GPU) ukazały się:

1. Wstęp teoretyczny.
2. Prototyp programu.
3. Pierwszy program GPU.
4. Program GPU v2.
5. Biblioteka CUBLAS. (Aktualnie czytany).

P.S.

Kod programu zamieszczam w załącznikach.

matmul_cublas.zip

[Qitr 41]

Jeśli w narodzie będzie wola ciągnięcia tematu to chętnie przedstawię jeszcze bibliotekę Thrust