Skalarni produkt dveh vektorjev¹

Poglejmo si sedaj skalarni produkt vektorjev. Začeli bomo s preprosto različico, da prikažemo osnovno idejo. Nato bomo analizirali učinkovitost preproste različice in jo razširili na različico, ki uporablja skupni pomnilnik.

Osnovna različica skalarnega produkta

The complete program for the naïve version of dot product is: Celotna koda programa za izračun skalarnega produkta

#include <stdio.h>
#include <stdlib.h>

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

#include "helper_cuda.h"

// CUDA kernel. Each thread takes care of one element of c
__global__ void dotProduct(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;
    }
}

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

    dim3 blockSize = 64;
    dim3 gridSize = n / blockSize.x;

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

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

    // Sum up the products
    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;
}

Izračun poteka v dveh korakih. V prvem najprej izračunamo produkte istoležnih elementov. Ščepec za izračun produkta:

// CUDA kernel. Each thread takes care of one element of c
__global__ void dotProduct(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;
    }
}

Rešitev je praktično enaka kot pri seštevanju dveh vektorjev, le da smo operacijo seštevanja nadomestili s produktom. V drugem koraku seštejemo vse produkte. Ko se ščepec zaključi, mora gostitelj prebrati vektor produktov iz naprave in sešteti vse njegove elemente v zanki:

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

Vse ostalo je enako kot pri primeru seštevanja vektorjev.

Zgornja koda se nahaja v mapi repozitorija delavnice skupaj z navodili za prevajanje in zagon 04-dot-product-naive.

Optimizacija skalarnega produkta: skupni pomnilnik in redukcija

Težava prejšnje implementacijo je, da gostitelj izračuna vsoto produktov parov. Pri zelo dolgih vektorjih lahko to vzame precej časa.

Zato bomo predstavili boljšo rešitev za seštevanje delnih produktov. Spoznali bomo, kako lahko niti uporabljajo skupni pomnilnik in sodelujejo pri seštevanju. Izračun skalarnega produkta bo sestavljen iz treh korakov:

• n niti bomo razdelili v b blokov niti s l nitmi na blok. Vsaka nit v bloku bo izračunala produkte istoležnih elementov vektrojev in jih seštela v svoj delni skalarni produkt v skupnem pomnilniku. Ker imamo l niti v bloku, bo na koncu tega koraka v deljenem pomnilniku shranjenih l delnih produktov.
• Vse niti bloka bodo nato sodelovale pri združevanju vsote l delnih skalarnih produktov. Ena izmed niti v bloku bo to vsoto shranila v globalni pomnilnik naprave. Za b blokov bomo imeli na koncu tega koraka b vsot v globalnem pomnilniku naprave.
• Sedaj bomo prenesli b vsot iz globalnega pomnilnika naprave na gostitelja in jih tam sešteli. Ker je število blokov b veliko manjše od velikosti vektorjev n, bo gostitelj moral prenesti bistveno manj podatkov iz naprave. Poleg tega bo gostitelj imel veliko manj dela pri končnem izračunu skalarnega produkta.

1. Korak: Izračun delnih skalarnih produktov

Spodnja koda prikazuje ščepec, ki zmnoži istoležne elemente, jih sešteje in delni skalarni produkt shrani v skupni pomnilnik (polje buf).

    // Get global thread ID
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    __shared__ float buf[BLOCK_SIZE];

    // Initialize shared buffer
    buf[threadIdx.x] = 0;

    // Make sure we do not go out of bounds
    while (tid < n)
    {
        buf[threadIdx.x] += a[tid] * b[tid];
        tid += gridDim.x * blockDim.x;
    }

    __syncthreads();

Tabela buf se hrani v skupnem pomnilniku in ima velikost BLOCK_SIZE. BLOCK_SIZE je konstanta in je enaka številu niti v bloku. S tem zagotovimo, da bo imela vsaka nit znotraj bloka prostor, kamor bo lahko zapisala delni skalarni produkt. Za rezervacijo skupnega pomnilnika preprosto uporabimo kvalifikator shared pri deklaraciji spremenljivke:

__shared__ float buf[BLOCK_SIZE];

Sledi inicializacija tabele buf z ničlami, saj se podatkovne strukture v skupnem pomnilniku ne inicializirajo samodejno. Tu, vsaka nit bloka inicializira svojo lokacijo v skupnem pomnilniku.

 buf[threadIdx.x] = 0;

Each thread can now multiply two corresponding elements and sums up the partial products into its entry in the buffer buf: Sedaj lahko vsaka nit zmnoži istoležne elemente in akumulira vsoto produktov na ustrezno mesto v tabeli buf.

    // Make sure we do not go out of bounds
    while (tid < n)
    {
        buf[threadIdx.x] += a[tid] * b[tid];
        tid += gridDim.x * blockDim.x;
    }

Pred začetkom drugega koraka združevanja delnih skalarnih produktov, moramo zagotoviti, da so vse niti znotraj bloka zaključile z izračuni do te točke. Na srečo obstaja funkcija:

 __syncthreads();

Ta klic zagotovi, da so vse niti znotraj bloka zaključile vse ukaze do tega trenutka in šele potem dovoli, da niti nadaljujejo z izvajanjem ukazov, ki sledijo. Temu ukazu rečemo prepreka (angl. barrier). Na splošno je prepreka sinhronizacijski konstrukt, ki ga uporabimo na določenem mestu v kodi in poskrbi, da nobena nit ne napreduje prek tega konstrukta, dokler vse niti (bloka) ne prispejo do njega. Ko vse niti dosežejo prepreko, se izvajanje nalednjih ukazov v kodi lahko nadaljuje. Če določene niti iz nekega razloga prepreke nikoli ne dosežejo, bodo preostale niti, ki so prepreko dosegle, preprosto obstale na prepreki. Prišlo bo do smrtnega objema (angl. dead-lock).

2. Korak: Redukcija delnih vsot

Zdaj, ko smo zagotovili, da je tabela buf napolnjena z delnimi skalarnimi produkti, lahko seštejemo vrednosti v njej. Proces, kjer vzamemo vhodno tabelo in izvedemo nekaj izračunov, ki proizvedejo manjše polje rezultatov, imenujemo redukcija. Naivni način izvedbe te redukcije bi bil, da uporabimo eno, ki se sprehodi preko tabele buf in računa tekočo vsoto. Časovna zahtevnost takega pristopa bi bila sorazmerna z dolžino polja. V bloku imamo na voljo lahko na stotine niti in lahko tako redukcijo izvedemo vzporedno dosežemo časovno zahtevnost, ki je sorazmerna z logaritmom dolžine polja.

Naslednja slika prikazuje postopek vzporedne redukcije:

Idea je, da vsaka nit vzame dve vrednosti iz tabele buf, ju sešteje in rezultat shrani nazaj v buf. Ker vsaka nit združi po dve vrednosti v eno, prvo stopnjo naredimo s polovico niti, s katerimi smo začeli. V naslednjem koraku isto stvar ponovimo s polovico niti iz prvega koraka. Nadaljujemo, dokler nimamo vsote vseh vrednosti tabele buf v prvem elementu.

Recimo, da imamo 16 elementov v polju buf in 16 niti v bloku. V prvem koraku vsaka od prvih osmih niti v bloku sešteje element z indeksom threadIdx.x in element z indeksom threadIdx.x + 8 ter rezultat zapiše v element z indeksom threadIdx.x. Spomnimo se, da je threadIdx.x indeks niti znotraj bloka. Preostalih osem niti iz bloka ne izvede ničesar. Nato vse niti počakajo pri prepreki, preden nadaljujejo z drugim korakom. S tem poskrbimo, da so vsa seštevanja zaključena, preden ponovno združujemo vrednosti.

V drugem koraku le prve štiri niti seštevajo elemente. Skupaj seštejemo elementa z indeksom threadIdx.x in indeksom threadIdx.x + 4, medtem ko druge niti ne izvedejo ničesar. Pred vstopom v tretji korak bodo vse niti znova počakale pri prepreki. Postopek nadaljujemo, dokler ne ostane le ena nit z indeksom 0, ki sešteje elementa z indeksom 0 in indeksom 1. Postopek redukcije, ki smo ga pravkar opisali, je prikazan v spodnjem izseku kode:

    // Reduction:
    int i = blockDim.x / 2;
    while (i != 0)
    {
        if (threadIdx.x < i)
        {
            // perform the addition:
            buf[threadIdx.x] += buf[threadIdx.x + i];
        }
        i = i / 2;
        // wait at the barrier:
        __syncthreads();
    }

Pomembno je opozoriti, da morajo vse niti v bloku klicati funkcijo __syncthreads(). Če se prepreka uporabi znotraj pogojnega stavka, moramo zagotoviti, da vse niti v bloku vstopijo v pogojni stavek. Na primer, spodnja koda prikazuje nepravilno uporabo prepreke, ker vse niti bloka ne vstopijo v pogojni stavek:

    if (threadIdx.x < i)
    {
        // perform the addition:
        buf[threadIdx.x] += buf[threadIdx.x + i];
        __syncthreads();
    }      

Niti z indeksom threadIdx.x, večjim ali enakim i, ne bodo nikoli izvedle klica __syncthreads(). Posledično bo del niti, za katere je bil pogoj izpolnjen ostal na prepreki v nedogled. To povzroči, da se ščepec "obesi". Tak ščepec bo povzročil, da se bo GPE se bo na neki točki prenehala odzivati in program bomo morali prisilno končati.

Po zaključku zanke while ima vsak blok niti na prvem mestu v tabeli buf seštevek vseh delnih skalarnih produtkov izračunanih znotraj bloka. Rezultat moramo sedaj zapisati v globalni pomnilnik, da ga bomo lahko kasneje prenesli na gostitelja. Ker moramo prenesti samo eno vrednost, uporabimo eno nit, ki vrednost zapiše na ustrezno mesto. To je določeno z indeksom bloka blockIdx.x:

    // only one thread in block writes the result:
    if (threadIdx.x == 0)
    {
        c[blockIdx.x] = buf[0];
    }

3. Korak: Izračun končnega rezultata na gostitelju

Zadnji korak zahteva, da na gostitelju seštejemo delne skalarne produkte, ki so jih izračunali posamezni bloki niti na napravi. Ker je polje c relativno majhno si to lahko brez težav privoščimo. Spodaj je prikazan koda celotnega programa za izračun skalarnega produkta s pomočjo skupnega pomnilnika in vzporedne redukcije.

#include <stdio.h>
#include <stdlib.h>

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

#include "helper_cuda.h"

#define BLOCK_SIZE 64

// CUDA kernel. Each thread takes care of one element of c, each thread block preforms reduction
__global__ void dotProduct(const float *a, const float *b, float *c, int n)
{
    // Get global thread ID
    int tid = blockIdx.x * blockDim.x + threadIdx.x;

    __shared__ float buf[BLOCK_SIZE];

    // Initialize shared buffer
    buf[threadIdx.x] = 0;

    // Make sure we do not go out of bounds
    while (tid < n)
    {
        buf[threadIdx.x] += a[tid] * b[tid];
        tid += gridDim.x * blockDim.x;
    }

    __syncthreads();

    // Reduction:
    int i = blockDim.x / 2;
    while (i != 0)
    {
        if (threadIdx.x < i)
        {
            // perform the addition:
            buf[threadIdx.x] += buf[threadIdx.x + i];
        }
        i = i / 2;
        // wait at the barrier:
        __syncthreads();
    }

    // only one thread in block writes the result:
    if (threadIdx.x == 0)
    {
        c[blockIdx.x] = buf[0];
    }
}

int main(int argc, char *argv[])
{
    // Size of vectors
    int n = 1024;
    dim3 blockSize = BLOCK_SIZE;
    dim3 gridSize = n / blockSize.x;

    size_t datasize = sizeof(float) * n;
    size_t partial_results_size = sizeof(float) * gridSize.x;

    // 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(partial_results_size);

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

    // 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));

    // Create CUDA events
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    // Record start event on the stream
    cudaEventRecord(start);

    // Execute the kernel
    dotProduct<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);

    // record stop event on the stream
    cudaEventRecord(stop);
    getLastCudaError("dotProduct() execution failed\n");

    // wait until the stop event completes
    cudaEventSynchronize(stop);    

    // Calculate the elapsed time between two events
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);
    printf("Kernel Execution time is: %0.3f milliseconds \n", milliseconds);

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

    // Sum up the products
    float result = 0.0;
    for (int i = 0; i < gridSize.x; 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));

    // Clean up the two events
    cudaEventDestroy(start);
    cudaEventDestroy(stop);

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

    return 0;
}

Zgornja koda se nahaja v mapi repozitorija delavnice skupaj z navodili za prevajanje in zagon 05-dot-product-reduction.

Merjenje časa izvajanja s pomočjo dogodkov CUDA

Napisali smo dve različici za izračun skalarnega produkta v okolju CUDA. Kako ugotovimo, katera različica je hitrejša? To lahko ugotovimo s pomočjo dogodkov CUDA v tokovih (angl. CUDA Streams).

Tok CUDA se nanaša na zaporedje operacij CUDA, ki se izvajajo na napravi v vrstnem redu, določenim s kodo gostitelja. Te operacije lahko vključujejo prenos podatkov med gostiteljem in napravo, zagon ščepca in večino drugih ukazov, ki jih izda gostitelj, vendar jih obdela GPU-naprava. Skozi to delavnico je tipičen vzorec v programiranju:

• Premakni vhodne podatke iz gostitelja na napravo.
• Izvedi ščepec na napravi.
• Premakni rezultat iz naprave nazaj na gostitelja.

Prenosi podatkov med gostiteljem in napravo s cudaMemcpy so sinhroni (ali blokirajoči). Sinhroni prenosi podatkov se ne začnejo, dokler niso bile zaključene vi prejšnji klici funkcij CUDA. Poleg tega se noben naslednji klic CUDA ne more začeti, dokler sinhroni prenos ni zaključen. Zagoni ščepcev pa so neblokirajoči. Ko ščepec zaženemo se izvajanje na CPE nadaljuje takoj in ne čaka na dokončanje ščepca. Prav blokirajoča narava prenosa podatkov preko funkcije cudaMemcpy zagotavlja, da se ščepec zaključi preden se izvrši prenos.

Dogodek CUDA je v bistvu časovni žig v toku CUDA, ki se zabeleži ob določenem trenutku med izvajanjem operacij v toku. Programski vmesnik CUDA nudi funkcije, ki vam omogočajo, da vstavite dogodke na katero koli mesto v toku in tudi preverite ali se je dogodek zgodil. Dogodek, ki ga vstavimo v nek tok, se bo zgodil šele, ko bodo vse predhodne operacije v istem toku zaključene. Uporaba dogodkov je relativno enostavna, saj zajemanje časovnega žiga vključuje le dva koraka: ustvarjanje dogodka in beleženje dogodka.

Dogodek dodamo v tok CUDA s pomočjo funkcije:

cudaError_t cudaEventRecord(cudaEvent_t event, cudaStream_t stream = 0);

Da se dogodek zgodi potem počakamo s pomočjo funkcije:

cudaError_t cudaEventSynchronize(cudaEvent_t event);

Zgornja fukcija bo blokirala izvajanje programa na gostitelju, dokler se vse predhodne operacije v nekem toku CUDA ne zaključijo.

S pomočjo sledeče funkcije lahko izmerimo čas med dvema dogodkoma v toku CUDA:

cudaError_t cudaEventElapsedTime(float* ms, cudaEvent_t start, cudaEvent_t stop);

Funkcija vrne pretečeni čas med dvema dogodkoma v milisekundah. This function returns the time elapsed between the events start and stop in milliseconds.

Spodnja koda prikazuje primer uporabe dogodkov za merjenje časa izvajanja ščepca dotProduct:

    // Create CUDA events
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    // Record start event on the stream
    cudaEventRecord(start);

    // Execute the kernel
    dotProduct<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);

    // Record stop event on the stream
    cudaEventRecord(stop);
    getLastCudaError("dot_product() execution failed\n");

    // Wait until the stop event completes
    cudaEventSynchronize(stop);    

    // Calculate the elapsed time between two events
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);
    printf("Kernel Execution time is: %0.3f milliseconds \n", milliseconds);

Zgornja koda se nahaja v mapi repozitorija delavnice skupaj z navodili za prevajanje in zagon 06-dot-product-reduction-profile.

---

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. ↩