Skoči na vsebino

Seštevanje vektorjev1

Tipično zaporedje operacij v programu, ki izkorišča GPE je naslednje:

  • Deklaracija in rezervacija pomnilnika na gostitelju in napravi.
  • Inicializacija podatkov na gostitelju.
  • Prenos podatkov iz gostitelja na napravo.
  • Izvedba enega ali več ščepcev.
  • Prenos rezultatov iz naprave na gostitelja.

Skozi primer programa, ki izračuna vsoto istoležnih elementov dveh vektorjev si poglejmo, kako izvedemo zgornje operacije s pomočjo vmesnika CUDA. Začnimo s preprosto funkcijo v programskem jeziku C, ki sešteje istoležne elemente dveh tabel (vektorjev) s po n elementi. Problem bomo rešili na dva načina. Najprej se bomo zadeve lotili precej naivno in se osredotočili samo na osnove programiranja za GPE in korake na gostitelju, ki so potrebni za zagon našega prvega programa na GPE. Nato bomo poskusili naš ščepec izboljšati.

Spodnja koda prikazuje implementacijo funkcije v C-ju za seštevanje dveh vektorjev, katerih elementi so števila s plavajočo vejico v enojni natančnosti:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// Add the elements of two arrays
void vecAdd(float *a,
               float *b,
               float *c,
               int n) {

    int tid = 0;

    while (tid < n) {
        c[tid] = a[tid] + b[tid];
        tid += 1;
    }
}

Argumenti funkcije so naslovi vseh treh vektorjev (a, b, c) in število elementov v vektorjih (n). Vektorja a in b sta seštevanca, medtem ko bomo v vektor c shranili rezultat. Seštevanje opravimo v zanki, ki se ponovi n-krat - toliko kot je dolžina vseh treh vektorjev.

Opazimo lahko, da so posamezni elementi vektorjev indeksirani z celoštevilskim indeksom tid. To ime smo izbrali namenoma. Razlog za takšno izbiro bomo kmalu razumeli.

Naivni pristop k seštevanju vektorjev na GPE

Sedaj bi radi funkcijo vectorAdd izvedli na GPE. Naučili smo se, da so GPE idealne za izvajanje enakih in neodvisnih operacij nad velikim številom podatkov. Naloga seštevanja vektorjev je očitno taka. Nad različnimi elementi vektorjev izvajamo isto operacijo, pri čemer so posamezna seštevanja med sabo neodvisna (opravimo jih lahko v poljubnem vrstnem redu ali pa celo vse naenkrat). Posledično je ta naloga zelo primerna za izvedbo na GPE.

Ščepec za seštevanje vektorjev

Da izvedemo seštevanjw vektorjev na GPE, jo moramo zapisati kot ščepec. Vsaka nit bo izvajala isto kodo. Osnovna ideja je nadomestiti iteracije zanke z instancami ščepca, ki se izvajajo nad vsakim elementom vektorja. Seštevanje dveh vektorjev dolžine n izvedemo s pomočjo n instanc ščepca, ki jih bo izvedlo n niti. Ščepec, ki ga bodo niti izvajale na GPE, je prikazano v spodnji kodi:

Delo, ki ga opravi funkcijo vectorAdd moramo torej razdeliti čim bolj enakomerno med niti. V našem primeru je to precej preprosta naloga, saj vsaka nit sešteva po dva istoležna elementa. Zato bomo na napravi pognali toliko niti, kolikor je elementov v vektorjih - v našem primeru je to n niti. Na podlagi indeksa niti določimo katera dva elementa bo neka nit seštela:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// CUDA kernel. Each thread takes care of one element of c
__global__ void vecAdd(const float *a,
                        const float *b,
                        float *c,
                        const int n)
{
    // Get global thread ID
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    // Make sure we do not go out of bounds
    if (tid < n)
        c[tid] = a[tid] + b[tid];
}

Vidimo, da je ščepec zelo izhodnišni funkciji v jeziku C, le da ima dodan ima kvalifikator __global__. Ta kvalifikator opozori prevajalnik, da naj funkcijo prevede za napravo ter, da jo bo klical gostitelj. Argumenti se ščepcu prenašajo enako kot običajni funkciji v C-ju.

Spremenljivke, definirane znotraj ščepca (npr. tid), se hranijo na napravi. Spremenljivke a, b in c morajo biti veljavni kazalci v pomnilniški prostor na napravi. To drži, ker bomo prostor za vektorje rezervirali v pomnilniku naprave, vektorja a in b pa bosta pred klicem ščepca prenesena iz pomnilnika gostitelja v pomnilnik naprave. Zadnjega argumenta, n, ni potrebno izrecno prenesti na napravo v kodi gostitelja. V jezikih C/C++ se argumenti funkcije privzeto prenašajo po vrednosti, za kar poskrbi izvajalno okolje CUDA.

Želimo, da vsaka nit sešteje en par istoležnih elementov vektorjev a in b ter shrani rezultat v vektor c. Indeks elementov, ki jih mora posamezna nit sešteti izračunamo na podlagi indeksa niti. Spomnimo se, nit je unikatno definirana s pomočjo indeksa bloka blockIdx in indeksa niti znoraj bloka threadIdx. V našem primeru, imamo oprava z vektorji, ki so enodimenzionalni, zato bomo tudi bloke in niti ob zagonu ščepca organizirali v eni dimenziji. Z izrazom:

int tid = blockIdx.x * blockDim.x + threadIdx.x;

izračunamo indeks elementov, ki ju bo posamezna nit seštela. Opazimo lahko, da ščepec ne vsebuje več zanke while. Razlog je preprosto to, da bo vsaka nit, ki izvaja ščepec izračunala vsoto dveh elementov, zagnali pa bomo vsaj toliko niti kot je elementov vektorja. Opazimo lahko, da preden izvedemo seštevanje preverimi, če je indeks tid znotraj dolžine vektorjev n. S tem poskrbimo za primere, ko število elementov vektroja ni deljivo z velikostjo bloka. V takem primeru moramo namreč zagnati več niti kot je dolžina vektorja.

Program na gostitelju

Aplikacija se izvaja na gostitelju in pošilja ščepce v izvajanje na prisotne naprave. Aplikacija na gostitelju mora izvesti naslednje korake:

  • Rezervira pomnilnik za vektorje v pomnilniku gostitelja.
  • Rezervira pomnilnik za vektorje v pomnilniku naprave.
  • Inicializira vektorja a in b v pomnilniku gostitelja.
  • Prenese vektorja a in b iz gostitelja na napravo.
  • Določi organizacijo niti.
  • Zažene ščepec.
  • Prenese vektor c iz naprave nazaj na gostitelja.

Spodnji je prikaza koda, ki se bo izvajala na gostitelju. V odstavkih, ki sledijo je pojasnjen vsak korak znotraj kode.

 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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
#include <stdio.h>
#include <stdlib.h>

#include <cuda.h>
#include <cuda_runtime.h>
#include "helper_cuda.h"

int main(int argc, char *argv[])
{
    // Size of vectors
    int n = 1024;
    size_t datasize = sizeof(float) * n;

    // Host input vectors
    float *h_a;
    float *h_b;
    // Host output vector
    float *h_c;

    // Device input vectors
    float *d_a;
    float *d_b;
    // Device output vector
    float *d_c;

    // Allocate memory for each vector on host
    h_a = (float *)malloc(datasize);
    h_b = (float *)malloc(datasize);
    h_c = (float *)malloc(datasize);

    // Allocate memory for each vector on GPU
    checkCudaErrors(cudaMalloc(&d_a, datasize));
    checkCudaErrors(cudaMalloc(&d_b, datasize));
    checkCudaErrors(cudaMalloc(&d_c, datasize));

    // Initialize vectors on host
    for (int i = 0; i < n; i++)
    {
        h_a[i] = 1.0;
        h_b[i] = 1.0;
    }

    // Copy host vectors to device
    checkCudaErrors(cudaMemcpy(d_a, h_a, datasize, cudaMemcpyHostToDevice));
    checkCudaErrors(cudaMemcpy(d_b, h_b, datasize, cudaMemcpyHostToDevice));

    // Define thread grid
    dim3 blockSize = 1024;
    dim3 gridSize = 1;

    // Execute the kernel
    vecAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);
    getLastCudaError("vecAdd() execution failed\n");

    // Copy array back to host
    checkCudaErrors(cudaMemcpy(h_c, d_c, datasize, cudaMemcpyDeviceToHost));

    // Check the result
    float result = 0.0;
    for (int i = 0; i < n; i++)
        result += h_c[i];
    printf("Result = %.2f \n", result);

    // Release device memory
    checkCudaErrors(cudaFree(d_a));
    checkCudaErrors(cudaFree(d_b));
    checkCudaErrors(cudaFree(d_c));

    // Release host memory
    free(h_a);
    free(h_b);
    free(h_c);

    return 0;
}

Znotraj funkcije main najprej definiramo vse kazalce za sklicevanje na vektorje v pomnilniku gostitelja in pomnilniku naprave

14
15
16
17
18
19
20
21
22
23
24
    // Host input vectors
    float *h_a;
    float *h_b;
    // Host output vector
    float *h_c;

    // Device input vectors
    float *d_a;
    float *d_b;
    // Device output vector
    float *d_c;

Rezerviramo prostor za vektorje na gostitelju, z argumentom datasize povemo želeno velikost pomnilnika v bajtih:

26
27
28
29
    // Allocate memory for each vector on host
    h_a = (float *)malloc(datasize);
    h_b = (float *)malloc(datasize);
    h_c = (float *)malloc(datasize);

Rezerviramo prostor za vektorje na napravi s pomočjo funckije cudaMalloc:

31
32
33
34
    // Allocate memory for each vector on GPU
    checkCudaErrors(cudaMalloc(&d_a, datasize));
    checkCudaErrors(cudaMalloc(&d_b, datasize));
    checkCudaErrors(cudaMalloc(&d_c, datasize));

Kazalci h_a, h_b in h_c kažejo v pomnilnik gostitelja, kazalci d_a, d_b in d_c kažejo v pomnilnik naprave. Pomnilnika gostitelja in naprave imata ločene pomnilniške prostore. Oba pa je mogoče upravljatiiz gostiteljske kode. Funkcija checkCudaErrors poskrbi za preverjanje napak pri izvajanju funkcij CUDA in izpiše ustrezno informacijo o napaki.

Gostitelj sedaj inicializira vektorja h_a in h_b:

36
37
38
39
40
41
    // Initialize vectors on host
    for (int i = 0; i < n; i++)
    {
        h_a[i] = 1.0;
        h_b[i] = 1.0;
    }

Po inicializaciji preprosto prenesemo podatke vektorjev h_a in h_b v ustrezne vektorje na napravi d_a in d_b z uporabo funkcije cudaMemcpy. Njen četrti argument določa smer kopiranja. V tem primeru uporabljamo cudaMemcpyHostToDevice, da določimo, da je prvi argument (cilj) kazalec na napravi, drugi argument (vir) pa kazalec na gostitelju:

43
44
45
    // Copy host vectors to device
    checkCudaErrors(cudaMemcpy(d_a, h_a, datasize, cudaMemcpyHostToDevice));
    checkCudaErrors(cudaMemcpy(d_b, h_b, datasize, cudaMemcpyHostToDevice));

Ena najpogostejših napak, ki jo naredijo tisti, ki se učijo programiranja v CUDA, je nepravilno sklicevanje do različnih pomnilniških prostorov. Na gostitelju ne smemo in ne moremo dostopati do podatkov preko kazalca na napravi, saj pomnilnik, rezerviran na GPE, ni neposredno naslovljiv s strani CPE. Novejše arhitekture GPE omogočajo dostop do pomnilnika naprave in gostitelja z uporabo enega samega kazalca preko t. i. koncepta enotnega pomnilnika (angl. Unified Memory), s katerim se pa na tej delavnici ne bomo ukvarjali.

Organizacija niti je ključni del programiranja v CUDA. Preden zaženemo ščepec, moramo nastaviti velikost mrežeo (tj. določiti dimenzijo bloka in dimenzijo mreže):

47
48
49
    // Define thread grid
    dim3 blockSize = 1024; // 1024 threads per block
    dim3 gridSize = 1; // one block

Končno, zaženemo ščepec:

51
52
    // Execute the kernel
    vecAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);
V našem primeru smo zagnali ščepec na 1024 nitih, ki so organizirane kot v enem bloku. V primeru, da so vektorji krajši od 1024, bo pogoj znotraj ščepca poskrbel, da niti ne bodo dostopale do pomnilniških lokacij izven vektorjev.

Ko se ščepec zaključi prenesemo rezultat nazaj v pomnilnik gostitelja: 

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
    // Copy array back to host
    checkCudaErrors(cudaMemcpy(h_c, d_c, datasize, cudaMemcpyDeviceToHost));

    // Check the result
    float result = 0.0;
    for (int i = 0; i < n; i++)
        result += h_c[i];
    printf("Result = %.2f \n", result);

    // Release device memory
    checkCudaErrors(cudaFree(d_a));
    checkCudaErrors(cudaFree(d_b));
    checkCudaErrors(cudaFree(d_c));

    // Release host memory
    free(h_a);
    free(h_b);
    free(h_c);

Zgornja koda je objavljena mapi repozitorija delavnice, skupaj z navodili za prevajanje in zagon 02-vec-add-short .

Seštevanje vektorjev poljubne dolžine

Med preučevanjem izvajalnega modela GPE smo opazili, da strojna oprema omejuje število blokov in število niti na blok, ki jih lahko zaženemo. Vsaka računska enota ima omejeno število registrov in omejeno količino skupnega pomnilnika. Kot smo že omenjali želimo imeti čim več aktivnih snopov med izvajanjem, saj s tem bolje izkoristimo GPE. Z visoko zasedenostjo računskih enot lahko pri izvajanju ukazov skrijemo zakasnitve. Spomnimo se, da je izvajanje drugih snopov, ko je en snop ustavljen, edini način, da zakrijemo zakasnitve in ohranimo zasedenost strojne opreme. Trenutno je naša koda omejena ne samo performančno ampak tudi funkcionalno. Seštevamo lahko le vektorje dolge največ 1024 elementov in pri tem izkoriščamo samo eno računsko enoto. Če hočemo narediti naš ščepec bolj splošen in zmogljiv bomo morali narediti nekaj sprememb in uporabiti več blokov niti pri računanju. 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
// CUDA kernel. Each thread takes care of one or more elements of c
__global__ void vecAdd( const float *a, 
                        const float *b, 
                        float *c, 
                        const int n)
{
    // Get global thread ID
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    // Make sure we do not go out of bounds
    while (tid < n){
        c[tid] = a[tid] + b[tid];
        tid += gridDim.x * blockDim.x; // Increase index by total number of threads in grid
    }
}

Ščepec, sedaj zgleda precej podoben funckiji iz katere smo izhajali na začetku. Ponovno uporabljamo zanko while, da se sprehajamo po vektorjih. Vendar, namesto, da indeks tid povečujemo za 1, ga povečujemo za skupno število niti, ki smo jih zagnali. Ko neka nit opravi seštevanje, ki ji je bilo dodeljeno, se premakne naprej po vektorjih in preveri, če slučajno obstaja še kak dodaten par elementov, ki ni bil dodeljen nobeni drugi niti. Ta skok naprej je enak za vse niti in je enostavno enak produktu dimenzij mreže in dimenzij bloka:

13
    tid += gridDim.x * blockDim.x;

Seveda, moramo tudi ustrezno spremeniti število organizacijo niti pri zagonu ščepca:

47
48
49
    // Define thread grid
    dim3 blockSize = 128; // 128 threads per block
    dim3 gridSize = n/blockSize; // compute the number of blocks needed

S tem enostavno razrešimo težavo kako lahko seštevamo vektorje poljubnih dolžin ne glede na to, kakšno organizacijo niti uporabimo. Mogoče, zaradi omejitev stroje opreme neke GPE, ne moremo zagnati dovolj niti, da bi vsaka dobila po en par elementov; mogoče želimo bolje izkoristiti računske enote določene GPE in zagnati večje število blokov z manj nitmi na blok. S prikazanim pristopom pokrijemo vse te primere. Zagotovo se sprašujete, zakaj ena nit preprosto ne sešteva 2 ali 4 sosednjih parov istoležnih elementov? Razlog leži v načinu dostopa niti do pomnilnika. Spomnimo se, da se dostop do pomnilnika izvaja v segmentih. Dve sosednji niti v snopu morata dostopati do dveh sosednjih pomnilniških besed, da zagotovimo združevanje pomnilniških dostopov.

Zgornja koda je objavljena v mapi repozitorija delavnice 03-vec-add-arb.

  1. © Patricio Bulić, Davor Sluga, Univerza v Ljubljani, Fakulteta za računalništvo in informatiko. Gradivo je objavljeno pod licenco Creative Commons Priznanje avtorstva-Nekomercialno-Deljenje pod enakimi pogoji 4.0 Mednarodna. ↩