Section author: Mia Doričić, Vedran Miletić

rocBLAS: ROCm Basic Linear Algebra Subprograms

Kod: https://github.com/ROCmSoftwarePlatform/rocBLAS

Dokumentacija: https://rocblas.readthedocs.io/

rocBLAS nudi za korištenje osnovne podprograme koji u domeni linearne algebre.

Postoji više tipova podataka koje rocBLAS koristi, a to su:

(za više informacija o definicijama i deklaracijama ovih tipova možete pogledati ovdje: https://rocblas.readthedocs.io/en/master/functions.html )

  • rocblas_int služi za određivanje koristi li se int32 (podatkovni model LP64) ili int64 (podatkovni model ILP64) (više informacija o 64-bitnim podatkovnim modelima)

  • rocblas_stride prolazi između matrica ili vektora u funkcijama tipa strided_batched

  • rocblas_half služi za predstavljanje 16-bitnog broja s pomičnim zarezom (tzv. polovična preciznost, engl. half-precision)

  • rocblas_bfloat16 služi za predstavljanje 16-bitnog Brainovog broja s pomičnim zarezom (više informacija o formatu s pomičnim zarezom bfloat16)

  • rocblas_float_complex služi za predstavljanje realnih i imaginarnih dijelova kompleksnog broja u zapisu s jednostrukom preciznošću

  • rocblas_double_complex služi za predstavljanje realnih i imaginarnih dijelova kompleksnog broja u zapisu s dvostrukom preciznošću

  • rocblas_handle je struktura koja sadrži kontekst biblioteke rocBLAS; mora se inicijalizirati putem funkcije rocblas_create_handle()

Relevantni konstantni rocBLAS parametri:

  • rocblas_operation služi za određivanje treba li matrica biti transponirana ili ne

  • rocblas_operation_none operira nad nepromijenjenom matricom

  • rocblas_operation_transpose operira sa transponiranom matricom

  • rocblas_operation_conjugate_transpose operira sa konjugatom transponirane matrice

Primjeri

Službeni primjer clients/samples/example_sgemm.cpp: https://github.com/ROCmSoftwarePlatform/rocBLAS/blob/develop/clients/samples/example_sgemm.cpp

Program prikazuje kako provesti množenje matrica u single float-precision formatu.

SGEMM stoji skraćeno za Single-precision floating-point General Matrix Multiply (pregled svih BLAS rutina.

Promotrimo funkciju main():

int main()
{
 rocblas_operation transa = rocblas_operation_none, transb = rocblas_operation_transpose;
 float             alpha = 1.1, beta = 0.9;

 rocblas_int m = DIM1, n = DIM2, k = DIM3;
 rocblas_int lda, ldb, ldc, size_a, size_b, size_c;
 int         a_stride_1, a_stride_2, b_stride_1, b_stride_2;
 rocblas_cout << "sgemm example" << std::endl;
 if(transa == rocblas_operation_none)
 {
     lda        = m;
     size_a     = k * lda;
     a_stride_1 = 1;
     a_stride_2 = lda;
     rocblas_cout << "N";
 }
 else
 {
     lda        = k;
     size_a     = m * lda;
     a_stride_1 = lda;
     a_stride_2 = 1;
     rocblas_cout << "T";
 }
 if(transb == rocblas_operation_none)
 {
     ldb        = k;
     size_b     = n * ldb;
     b_stride_1 = 1;
     b_stride_2 = ldb;
     rocblas_cout << "N: ";
 }
 else
 {
     ldb        = n;
     size_b     = k * ldb;
     b_stride_1 = ldb;
     b_stride_2 = 1;
     rocblas_cout << "T: ";
 }
 ldc    = m;
 size_c = n * ldc;

 // Naming: da is in GPU (device) memory. ha is in CPU (host) memory
 std::vector<float> ha(size_a);
 std::vector<float> hb(size_b);
 std::vector<float> hc(size_c);
 std::vector<float> hc_gold(size_c);

 // initial data on host
 srand(1);
 for(int i = 0; i < size_a; ++i)
 {
     ha[i] = rand() % 17;
 }
 for(int i = 0; i < size_b; ++i)
 {
     hb[i] = rand() % 17;
 }
 for(int i = 0; i < size_c; ++i)
 {
     hc[i] = rand() % 17;
 }
 hc_gold = hc;

 // allocate memory on device
 float *da, *db, *dc;
 CHECK_HIP_ERROR(hipMalloc(&da, size_a * sizeof(float)));
 CHECK_HIP_ERROR(hipMalloc(&db, size_b * sizeof(float)));
 CHECK_HIP_ERROR(hipMalloc(&dc, size_c * sizeof(float)));

 // copy matrices from host to device
 CHECK_HIP_ERROR(hipMemcpy(da, ha.data(), sizeof(float) * size_a, hipMemcpyHostToDevice));
 CHECK_HIP_ERROR(hipMemcpy(db, hb.data(), sizeof(float) * size_b, hipMemcpyHostToDevice));
 CHECK_HIP_ERROR(hipMemcpy(dc, hc.data(), sizeof(float) * size_c, hipMemcpyHostToDevice));

 rocblas_handle handle;
 CHECK_ROCBLAS_ERROR(rocblas_create_handle(&handle));

 CHECK_ROCBLAS_ERROR(
     rocblas_sgemm(handle, transa, transb, m, n, k, &alpha, da, lda, db, ldb, &beta, dc, ldc));

 // copy output from device to CPU
 CHECK_HIP_ERROR(hipMemcpy(hc.data(), dc, sizeof(float) * size_c, hipMemcpyDeviceToHost));

 rocblas_cout << "m, n, k, lda, ldb, ldc = " << m << ", " << n << ", " << k << ", " << lda
              << ", " << ldb << ", " << ldc << std::endl;

 float max_relative_error = std::numeric_limits<float>::min();

 // calculate golden or correct result
 mat_mat_mult<float>(alpha,
                     beta,
                     m,
                     n,
                     k,
                     ha.data(),
                     a_stride_1,
                     a_stride_2,
                     hb.data(),
                     b_stride_1,
                     b_stride_2,
                     hc_gold.data(),
                     1,
                     ldc);

 for(int i = 0; i < size_c; i++)
 {
     float relative_error = (hc_gold[i] - hc[i]) / hc_gold[i];
     relative_error       = relative_error > 0 ? relative_error : -relative_error;
     max_relative_error
         = relative_error < max_relative_error ? max_relative_error : relative_error;
 }
 float eps       = std::numeric_limits<float>::epsilon();
 float tolerance = 10;
 if(max_relative_error != max_relative_error || max_relative_error > eps * tolerance)
 {
     rocblas_cout << "FAIL: max_relative_error = " << max_relative_error << std::endl;
 }
 else
 {
     rocblas_cout << "PASS: max_relative_error = " << max_relative_error << std::endl;
 }

 CHECK_HIP_ERROR(hipFree(da));
 CHECK_HIP_ERROR(hipFree(db));
 CHECK_HIP_ERROR(hipFree(dc));
 CHECK_ROCBLAS_ERROR(rocblas_destroy_handle(handle));
 return EXIT_SUCCESS;
}

Program odmah počinje sa deklaracijom dvije varijable; transa i transb, te su obje tipa rocblas_operation, iz razloga što su njihove vrijednosti izjednačene sa rocBLAS parametrima rocblas_operation_none i rocblas_operation_transpose. (Za više informacija o ovim parametrima pogledajte početak ovog teksta.)

Zatim se postavljaju varijable alpha i beta tipa float na određene vrijednosti:

rocblas_operation transa = rocblas_operation_none, transb = rocblas_operation_transpose;
float             alpha = 1.1, beta = 0.9;

Također se definiraju tri dimenzije koje će se koristiti, u obliku varijabli m, n i k tipa rocblas_int. Isto tako, definiraju se i vodeće dimenzije za A, B i C te njihove veličine, također tipa rocblas_int. (Za više informacija o ovom tipu pogledajte na početak ovog teksta).

rocblas_int m = DIM1, n = DIM2, k = DIM3;
rocblas_int lda, ldb, ldc, size_a, size_b, size_c;

Nakon toga slijedi definicija int varijabli stride za A i B, koje označavaju broj lokacija u memoriji između elemenata nekog polja. U ovom slučaju imamo:

int  a_stride_1, a_stride_2, b_stride_1, b_stride_2;

Ispisuje se putem rocblas_cout-a da je ovo primjer SGEMM raučunanja, te krećemo na if provjeru, gdje se kao glavni uvjet transa uspoređuje s rocblas_operation_none, te provjerava vrijedi li ta usporedba.

Ako vrijedi, tada:

  • vodeća dimenzija za A postaje m (dakle DIM1)

  • veličina za A postaje jednaka umnošku k i lda

  • stride_1 od A se postavlja na 1

  • stride_2 od A se postavlja na vrijednost lda

  • putem rocblas_cout program vraća slovo N, potvrđujući zadani uvjet.

rocblas_cout << "sgemm example" << std::endl;
if(transa == rocblas_operation_none)
{
    lda        = m;
    size_a     = k * lda;
    a_stride_1 = 1;
    a_stride_2 = lda;
    rocblas_cout << "N";
}

Ako je slučaj imalo drugačiji od zadanog uvjeta u if, tada:

  • vodeća dimenzija za A postaje k (dakle DIM3)

  • veličina za A postaje jednaka umnošku m i lda

  • stride_1 od A se postavlja na vrijednost lda

  • stride_2 od A se postavlja na 1

  • putem rocblas_cout-a program vraća slovo T, potvrđujući razliku od zadanog uvjeta.

else
{
   lda        = k;
   size_a     = m * lda;
   a_stride_1 = lda;
   a_stride_2 = 1;
   rocblas_cout << "T";
}

Ako pogledamo sljedeći if, primjetiti ćemo da je princip jednak kao i prijašnji, samo što se koriste varijable vezane za B. Ovdje kao zadani uvjet uspoređujemo transb sa rocblas_operation_none:

if(transb == rocblas_operation_none)
{
    ldb        = k;
    size_b     = n * ldb;
    b_stride_1 = 1;
    b_stride_2 = ldb;
    rocblas_cout << "N: ";
}
else
{
    ldb        = n;
    size_b     = k * ldb;
    b_stride_1 = ldb;
    b_stride_2 = 1;
    rocblas_cout << "T: ";
}

Slijedi postavljanje vodeće dimezije za C, te ju postavljamo na m (DIM1). Također, njenu veličinu računamo kao umnožak DIM2 i ldc.

ldc    = m;
size_c = n * ldc;

Inicijaliziramo vektore, gdje za svaki vektor pišemo slovo ‘h’ ispred, koje označava da je u memoriji procesora, host memory.

std::vector<float> ha(size_a);
std::vector<float> hb(size_b);
std::vector<float> hc(size_c);
std::vector<float> hc_gold(size_c);

Koristeći srand, kroz tri for petlje postavljaju se početne vrijednosti na svaki element memorije procesora, za A, B i C. Nakon svih petlji varijabla hc_gold izjednačena je sa dodjeljenim vrijednostima od hc.

srand(1);
for(int i = 0; i < size_a; ++i)
{
    ha[i] = rand() % 17;
}
for(int i = 0; i < size_b; ++i)
{
    hb[i] = rand() % 17;
}
for(int i = 0; i < size_c; ++i)
{
    hc[i] = rand() % 17;
}
hc_gold = hc;

Alociramo memoriju na uređaju putem funkcije hipMalloc i provjeravamo uspješnost korištenjem pomoćne funkcije CHECK_HIP_ERROR:

float *da, *db, *dc;
CHECK_HIP_ERROR(hipMalloc(&da, size_a * sizeof(float)));
CHECK_HIP_ERROR(hipMalloc(&db, size_b * sizeof(float)));
CHECK_HIP_ERROR(hipMalloc(&dc, size_c * sizeof(float)));

Zatim se matrice sa domaćina (CPU) kopiraju na uređaj, također pritom provjeravajući hoće li se pojaviti greška:

CHECK_HIP_ERROR(hipMemcpy(da, ha.data(), sizeof(float) * size_a, hipMemcpyHostToDevice));
CHECK_HIP_ERROR(hipMemcpy(db, hb.data(), sizeof(float) * size_b, hipMemcpyHostToDevice));
CHECK_HIP_ERROR(hipMemcpy(dc, hc.data(), sizeof(float) * size_c, hipMemcpyHostToDevice));

Stvara se rocBLAS handle, odnosno poveznica do rocBLAS biblioteke, koju stvaramo koristeći provjeru je li se putem pojavila greška. Nakon što smo stvorili handle, pozivamo funkciju rocblas_sgemm sa svim parametrima koje smo vidjeli definiranje u main() funkciji:

CHECK_ROCBLAS_ERROR(rocblas_create_handle(&handle));

CHECK_ROCBLAS_ERROR(
    rocblas_sgemm(handle, transa, transb, m, n, k, &alpha, da, lda, db, ldb, &beta, dc, ldc));

Rezultat ove funkcije kopiramo sa uređaja (GPU) na CPU uz HIP_ERROR provjeru, te ispisujemo dobivene vrijednosti:

CHECK_HIP_ERROR(hipMemcpy(hc.data(), dc, sizeof(float) * size_c, hipMemcpyDeviceToHost));

rocblas_cout << "m, n, k, lda, ldb, ldc = " << m << ", " << n << ", " << k << ", " << lda
             << ", " << ldb << ", " << ldc << std::endl;

Slijedi inicijalizacija varijable max_relative_error tipa float koja će nam trebati u sljedećim koracima, međutim, postavljeno joj je ograničenje:

float max_relative_error = std::numeric_limits<float>::min();

Funkcija min() ograničava ovu varijablu na način da postavlja njen minimum na najnižu moguću pozitivnu vrijednost, te spriječava tu vrijednost da ne ode u negativ (za više informacija ovoj funkciji proučite std::numeric_limits::min na cppreference.com).

Došli smo do poziva funkcije mat_mat_multi(...) (Matrix Matrix Multiply), koja je definirana na početku koda, a izgleda ovako:

void mat_mat_mult(T   alpha,
              T   beta,
              int M,
              int N,
              int K,
              T*  A,
              int As1,
              int As2,
              T*  B,
              int Bs1,
              int Bs2,
              T*  C,
              int Cs1,
              int Cs2)
{
for(int i1 = 0; i1 < M; i1++)
{
    for(int i2 = 0; i2 < N; i2++)
    {
        T t = 0.0;
        for(int i3 = 0; i3 < K; i3++)
        {
            t += A[i1 * As1 + i3 * As2] * B[i3 * Bs1 + i2 * Bs2];
        }
        C[i1 * Cs1 + i2 * Cs2] = beta * C[i1 * Cs1 + i2 * Cs2] + alpha * t;
    }
}
}

U njenoj definiciji možemo vidjeti inicijalizaciju svih potrebnih varijabli za dimenzije, veličine i sl.

Slijede tri for petlje koje uvjetuju jedna drugu; dok svaki brojač za svaku dimenziju vrti petlju do maksimalne vrijednosti svake dimezije, provodi se uvećavanje t varijable, čiji je iznos definiran umnoškom A i B matrica, čije su veličine zbroj umnožaka vrijednosti brojača i veličine trenutnog elementa. Isti princip koristi se za matricu C koja se izračunava na kraju druge for petlje.

for(int i1 = 0; i1 < M; i1++)
 {
     for(int i2 = 0; i2 < N; i2++)
     {
         T t = 0.0;
         for(int i3 = 0; i3 < K; i3++)
         {
             t += A[i1 * As1 + i3 * As2] * B[i3 * Bs1 + i2 * Bs2];
         }
         C[i1 * Cs1 + i2 * Cs2] = beta * C[i1 * Cs1 + i2 * Cs2] + alpha * t;
     }
 }

Nakon poziva ove funkcije u main()-u, program nastavlja sa petljom for u kojoj se postavlja vrijednost prije-spomenutog relative_error-a, nakon kojeg se izračunava i max_relative_error putem usporedbi veličina ove dvije varijable (za više informacija o usporedbi upotrebom ? znaka, odnosno “conditional operator”-a, pogledajte ovdje: https://en.cppreference.com/w/cpp/language/operator_other). Zatim se i koristi vrijednost epsilon-a, koji već postoji u biblioteci numeric_limits (za više informacija o ovoj funkciji pogledajte std::numeric_limits::epsilon na cppreference.com). Slijedi ispis, ovisno o ishodu operacije, te oslobađanje prethodno alocirane memorije.

for(int i = 0; i < size_c; i++)
{
    float relative_error = (hc_gold[i] - hc[i]) / hc_gold[i];
    relative_error       = relative_error > 0 ? relative_error : -relative_error;
    max_relative_error
        = relative_error < max_relative_error ? max_relative_error : relative_error;
}
float eps       = std::numeric_limits<float>::epsilon();
float tolerance = 10;
if(max_relative_error != max_relative_error || max_relative_error > eps * tolerance)
{
    rocblas_cout << "FAIL: max_relative_error = " << max_relative_error << std::endl;
}
else
{
    rocblas_cout << "PASS: max_relative_error = " << max_relative_error << std::endl;
}

CHECK_HIP_ERROR(hipFree(da));
CHECK_HIP_ERROR(hipFree(db));
CHECK_HIP_ERROR(hipFree(dc));
CHECK_ROCBLAS_ERROR(rocblas_destroy_handle(handle));
return EXIT_SUCCESS;

Kraj programa.

Pogledajmo drugi primjer upotrebe rocBLAS funkcija na drugom programu.

Službeni primjer clients/samples/example_sscal.cpp: https://github.com/ROCmSoftwarePlatform/rocBLAS/blob/develop/clients/samples/example_sscal.cpp

U ovom programu vidimo primjer skaliranja svakog elementa vektora x sa skalarom alpha.

SSCAL stoji skraćeno za Single-precision floating-point Scalar.

Program izgleda ovako:

int main()
{
 rocblas_int    N      = 10240;
 rocblas_status status = rocblas_status_success;
 float          alpha  = 10.0;

 std::vector<float> hx(N);
 std::vector<float> hz(N);
 float*             dx;

 double gpu_time_used;

 rocblas_handle handle;
 rocblas_create_handle(&handle);

 hipMalloc(&dx, N * sizeof(float));

 srand(1);
 rocblas_init<float>(hx, 1, N, 1);

 hz = hx;

 hipMemcpy(dx, hx.data(), sizeof(float) * N, hipMemcpyHostToDevice);

 printf("N        rocblas(us)     \n");

 gpu_time_used = get_time_us_sync_device(); // in microseconds


 status = rocblas_sscal(handle, N, &alpha, dx, 1);
 if(status != rocblas_status_success)
 {
     return status;
 }

 gpu_time_used = get_time_us_sync_device() - gpu_time_used;

 hipMemcpy(hx.data(), dx, sizeof(float) * N, hipMemcpyDeviceToHost);

 bool error_in_element = false;
 for(rocblas_int i = 0; i < N; i++)
 {
     if(hz[i] * alpha != hx[i])
     {
         printf("error in element %d: CPU=%f, GPU=%f ", i, hz[i] * alpha, hx[i]);
         error_in_element = true;
         break;
     }
 }

 printf("%d    %8.2f\n", (int)N, gpu_time_used);

 if(error_in_element)
 {
     printf("SSCAL TEST FAILS\n");
 }
 else
 {
     printf("SSCAL TEST PASSES\n");
 }

 hipFree(dx);
 rocblas_destroy_handle(handle);
 return rocblas_status_success;
}

Krenimo promatrati funkciju main().

Na početku se definiraju potrebne varijable te im se dodjeljuju vrijednosti. Pojavljuju se tipovi rocblas_int za varijablu N (za više informacija o ovom tipu pogledajte na početak teksta), te rocblas_status za varijablu status koja je izjednačena sa rocBLAS parametrom rocblas_status_success (za više informacija o ovom parametru pogledajte početak teksta). Skalar alpha ima vrijednost postavljenu na 10.0.

rocblas_int    N      = 10240;
rocblas_status status = rocblas_status_success;
float          alpha  = 10.0;

Slijedi inicijalizacija dva vektora hx i hz (h označava memoriju domaćina (engl. host), ondosno memoriju procesora), te pokazivača tipa float naziva dx (memorija uređaja (engl. device), odnosno memoriju grafičkog procesora). Također, stvara se varijabla gpu_time_used koja se koristi u kasnijim koracima za dohvaćanje vremena izvođenja operacija.

std::vector<float> hx(N);
std::vector<float> hz(N);
float*             dx;

double gpu_time_used;

Stvara se rocBLAS handle, odnosno poveznica do rocBLAS biblioteke, koju stvaramo koristeći provjeru je li se putem pojavila greška. Nakon što smo stvorili handle, alocira se memorija na uređaju (GPU).

rocblas_handle handle;
rocblas_create_handle(&handle);

hipMalloc(&dx, N * sizeof(float));

Koristeći srand, program čita početne podatke na procesoru. Također, kopira se vrijednost vektora hz u hx, koju zatim kopiramo sa CPU na GPU, te ispisujemo dobivenu vrijednost.

srand(1);
rocblas_init<float>(hx, 1, N, 1);

hz = hx;

hipMemcpy(dx, hx.data(), sizeof(float) * N, hipMemcpyHostToDevice);

Dohvaćamo vrijeme koje je bilo potrebno za operaciju putem funkcije get_time_us_sync_device() u mikrosekundama.

gpu_time_used = get_time_us_sync_device();

Varijabla status bit će jednaka ishodu funkcije rocblas_sscal() u kojoj se koriste prije-definirani handle, N, adresa od alpha i dx. Pokreće se if uvjetovanje koji vraća vrijednost varijable status ako i samo ako status ne pokazuje uspjeh putem parametra rocblas_status_success.

status = rocblas_sscal(handle, N, &alpha, dx, 1);
if(status != rocblas_status_success)
{
    return status;
}

Računa se vrijeme potrebno za izvođenje operacije do ovog koraka, na način da se oduzima vrijednost koju smo računali prethodno sa istom varijablom. Zatim, kopira se ispis podataka sa uređaja (GPU) na procesor.

gpu_time_used = get_time_us_sync_device() - gpu_time_used;

hipMemcpy(hx.data(), dx, sizeof(float) * N, hipMemcpyDeviceToHost);

Nakon svega, provjerava se je li uspješno provedena računica funkcije rocblas_scal. Na početku je vrijednost varijable koja označava moguću grešku u operaciji (tipa bool) na false. Ako dodđe do greške, ta vrijednost se mijenja u true. Slijedi ispis vremena potrebnog za izvršavanje operacije.

bool error_in_element = false;
for(rocblas_int i = 0; i < N; i++)
{
    if(hz[i] * alpha != hx[i])
    {
        printf("error in element %d: CPU=%f, GPU=%f ", i, hz[i] * alpha, hx[i]);
        error_in_element = true;
        break;
    }
}

printf("%d    %8.2f\n", (int)N, gpu_time_used);

Za kraj primjera, ispisuje se da, u slučaju da se kroz prijašnju for petlju pojavila greška, je izračun bio neuspješan. U suprotnom, ispisuje se da je test bio uspješan. Oslobađa se prethodno alocirana memorija, te se handle uništava.

if(error_in_element)
{
    printf("SSCAL TEST FAILS\n");
}
else
{
    printf("SSCAL TEST PASSES\n");
}

hipFree(dx);
rocblas_destroy_handle(handle);

return rocblas_status_success;

Kraj programa.