Množenje matrik¹

Sedaj se lotimo problema, ki nam bo pomagal razumeti 2-dimenzionalni prostor NDRange ter uporabo lokalnega pomnilnika za pohitritev izvajanj ščepcev na napravah GPE. V ta namen bomo sprogramirali množenje matrik. Brez izgube splošnosti, se bomo omejili na kvadratne matrike velikosti NxN. Množenje kvadratnih matrik A in B prikazuje spodnja slika:

Element matrike C[i,j] dobimo tako, da skalarno zmnožimo i-to vrstico matrike A ter j-ti stolpec matrike B. Torej, za vsak element matrike C moramo izračunati en skalarni prodiukt med vrstico iz A ter stolpcem iz B. Zgornja slika prikazuje, katere vrstice in stolpce je medseboj treba skalarno zmnožiti da dobimo elemente C[1,2], C[5,2] in C[5,7]. Na primer, za C[5,7] moramo skalarno zmnožiti 5. vrstico iz matrike A ter 7. stolpec iz matrike B.

Naivno množenje matrik

S preprosto verzijo množenja matrik, bomo najprej spoznali uporabo dvo-dimenzionalnih prostorov NDRange.

Ščepec

Tokrat se še ne bomo ukvarjali z učinkovitostjo ščepca. Pripravili ga bomo tako, da bo vsaka nit, ki ga bo izvajala, izračunala le en element matrike produkta.

__kernel void matmul_naive(
                __global float* matA,
                __global float* matB,
                __global float* matC,
                int N){


    int i = get_global_id(1);    // row (Y)
    int j = get_global_id(0);    // column (X)

    float dp;

    dp = 0.0;
    for (int k = 0; k < N; k++){
        dp += matA[i*N + k] * matB[k*N + j];
    }

    matC[i*N + j] = dp;
}

Ker so elementi matrik organizirani v 2-dimenzionalno polje, je najbolj naravno, da tudi niti organiziramo v dvo-dimenzionalni prostor NDRange. Vsaka nit bo sedaj določena s parom globalnik indeksov (i,j). V zgornji kodi zato vsaka nit najprej pridobi svoj globalni indeks s klicem funkcije get_global_id (vrstici 8 in 9).

    int i = get_global_id(1);    // row (Y)
    int j = get_global_id(0);    // column (X)

Nato vsaka nit izračuna skalarni produkt med i-to vrstico v matriki A in j-tim stolpcem v matriki B in rezultat shrani v element C[i,j] matrike produkta (vrstice 13 do 18).

Program na gostitelju

Program na gostitelju mora rezervirati prostor za vse tri matrike na gostitelju in inicializirati matriki A in B ter rezervirati prostor v globalnem pomnilniku naprave za tri matrike velikosti NxN (vrstice 20 do 34)

    cl_float* matA_h;
    cl_float* matB_h;
    cl_float* matC_h;

    // Size of data:
    size_t datasize = sizeof(cl_float) * iNumElements * iNumElements;

    // Allocate host arrays
    matA_h = (void *)malloc(datasize);
    matB_h = (void *)malloc(datasize);
    matC_h = (void *)malloc(datasize);
    // init arrays:
    for (int i = 0; i<iNumElements; i++ ) {
        for (int j = 0; j<iNumElements; j++ ) {
            matA_h[i*iNumElements + j] = 1.0;
            matB_h[i*iNumElements + j] = 1.0;
        }
    }

    cl_mem matA_d; // Input array on the device
    cl_mem matB_d; // Input array on the device
    cl_mem matC_d; // Output array on the device
    s
    // Use clCreateBuffer() to create a buffer object (d_A)
    // that will contain the data from the host array A
    matA_d = clCreateBuffer(
                             context,
                             CL_MEM_READ_ONLY,
                             datasize,
                             NULL,
                             &status);
    clerr_chk(status);

    // do the same for matB_d  an d matC_d ...

Poleg tega mora program na gostitelju nastaviti globalno velikost prostora NDRange ter število niti v delovni skupini.

    const size_t szLocalWorkSize[2] = {16, 16};
    const size_t szGlobalWorkSize[2] = {N, N};

V zgornjem primeru smo določili, da bomo zagnali NxN niti v dvo-dimenzionalnem prostoru NDRange ter da bo vsaka delovna skupina vsebovala 16x16 niti.

Ščepec na napravi GPE zaženemo s funkcijo clEnqueueNDRangeKernel, kjer s tretjim argumentom določimo, da je prostor NDRange dvo-dimenzionalen.

    status = clEnqueueNDRangeKernel(
                                   cmdQueue,
                                   ckKernel,
                                   2,
                                   NULL,
                                   szGlobalWorkSize,
                                   szLocalWorkSize,
                                   0,
                                   NULL,
                                   NULL);

Množenje matrik po ploščicah s pomočjo lokalnega pomnilnika

Prejšnja koda za množenje matrik ima kar nekaj hudih pomankljivosti. Vsaka nit namreč iz globalnega pomnilnika bere 2xN elementov: vrstico iz matrike A in stolpec iz matrike B. Pri pozornem opazovanju vidimo, da vse niti, ki računajo element matrike C v isti vrstici dostopajo do iste vrstice matrike A, vendar vsaka do različnega stolpca matrike B. Najhuje pri tem je, da pri dostopu do stolpcev v matirki B ne upoštevamo usklajenga dostopa do globalnega pomnilnika. Sosednje niti namreč ne dostopajo do sosednjih elementov matrike B, kar moćno upočasni izvajanje ščepca, saj sta dva sosednja elementa v stolpcu pri dolgih vrsticah lahko narazen za več kot 128 bajtov ter padeta v različne segmente. V takem primeru je za vsak element stolpca iz matrike B poteben ločen dostop do pomnilnika!

Po drugi strani pa je znano, da množenje matrik lahko izvedemo po ploščicah - množimo manjše podmatrike (ploščice), kot je prikazano na spodnji sliki.

Isti rezultat bomo torej dobili, če bomo elemente podmatrike C računali postopoma z množenjem podmatrik (ploščic) matrike A in matrike B. Sledeč iz zgornje slike je ideja sledeča: V lokalni pomnilnik najprej naložimo svetlozeleno ploščico matrike A in svetlomodro ploščico matrike B ter ju medseboj zmožimo. Tako dobimo delno izračunano rdečo ploščico matrike C. V naslednjem koraku iz matrike A naložimo temnozeleno ploščico ter iz matrike B temnomodro ploščico, ju medseboj zmnožimo ter rezultat prištejemo k rdeči ploščici iz matrike C. Ker sta v nekem koraku obe ploščici iz matrik A in B v lokalnem spominu, nimamo težav z usklajenim dostopom do globalnega pomnilnika. Povrh pa je dostop do lokalnega pomnilnika 100-krat hitrejši kot dostop do globalnega pomnilnika!

Ščepec, ki množi matriki po ploščicah predstavlja spodnja koda.

#define TW 16

__kernel void matmul_tiled( __global float* matA,
                            __global float* matB,
                            __global float* matC,
                            int N){

    // Local memory to fit the tiles
    __local float tileA[TW][TW];
    __local float tileB[TW][TW];

    // global thread index
    int rowGID = get_global_id(1);    // row in NDRange
    int colGID = get_global_id(0);    // column in NDRange

    // local thread index
    int rowLID = get_local_id(1);    // row in tile
    int colLID = get_local_id(0);    // column in tile

    float dp = 0.0;

    // Iterate through all the tiles necessary to calculate the products:
    for (int iTile = 0; iTile < N/TW; iTile++){
        // Collaborative load of tiles:
        // - Load a tile of matrix A into local memory:
        tileA[rowLID][colLID] = matA[rowGID*N + (colLID + iTile*TW)];
        // - Load a tile of matrix B into local memory:
        tileB[rowLID][colLID] = matB[(rowLID + iTile*TW)*N+colGID];

        // Synchronise to make sure the tiles are loaded
        barrier(CLK_LOCAL_MEM_FENCE);

        for(int k = 0; k < TW; k++) {
            dp += tileA[rowLID][k] * tileB[k][colLID];
        }
        // Wait for all WI to finish before loading next tile:
        barrier(CLK_LOCAL_MEM_FENCE);
    }

    // save the dot product into matrix C in global memory:
    matC[rowGID*N + colGID] = dp;
}

Vsaka nit se v zanki for sprehodi po vseh ploščicah, potrebnih za izračun njenega elementa matrike C (vrstice 23 do 38).

for (int iTile = 0; iTile < N/TW; iTile++){
    ...
}

V vsaki iteraciji zanke for niti iz iste delovne skupine najprej prenesejo ploščici iz matrike A in matrike B v lokalni pomnilnik ter se počakajo ob prepreki (vrstice 24 do 31).

    // Collaborative load of tiles:
        // Load a tile of matrix A into local memory:
    tileA[rowLID][colLID] = matA[rowGID*N + (colLID + iTile*TW)];
        // Load a tile of matrix B into local memory:
    tileB[rowLID][colLID] = matB[(rowLID + iTile*TW)*N+colGID];

    // Synchronise to make sure the tiles are loaded
        barrier(CLK_LOCAL_MEM_FENCE);
}

Tukaj se moramo za hipec ustaviti ter razložiti indekse, ki jih uporabimo pri branu ploščic iz globalnega pomnilnika naprave. Predpostavimo, da so ploščice dimenzije TWxTW, pri čemer je TW=4, ter da je N=8. Predpostavimo, da nit računa element C[5,2]. Ta element sodi v ploščico z indeksom (1,0). Nit, ki je zadožena za izračun elementa C[5,2] ima globalni indeks (5,2) in lokalni indeks (1,2). Ta nit mora sedaj dostopati do ploščic (1,0) in (1,1) iz matrike A in ploščic (0,0) and (1,1) iz matrike B. Nit (5,2) bo zadolžena za nalaganje naslednjih elementov v lokalni pomnilnik:

1. A(5,2) in B(1,2) v prvi iteraciji ter
2. A(5,6) in B(5,2) v drugi iteraciji.

Spodnja slika prikazuje relacijo med lokalnim indeksom posameznega elementa iz ploščice, ter lokalnim in globalnim indeksom niti, ki ga bere.

Koda na gostitelju bo ostala enaka kot pri naivnem množenju matrik, le naložiti moramo ščepec matmul_tiled.

Celotno kodo iz tega poglavja najdete v mapi 08-matrix-mul tukaj.

---

---

---

1. © Patricio Bulić, 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. ↩