funkcionálně.cz

Přední český blog o funkcionálním programování, kde se o funkcionálním programování nepíše
««« »»»

Hořící křemík & násobení matic

15. 6. 2017 — k47

Mám ten­denci jít vždycky příliš hlu­boko.

Na­pří­klad moje ne­dávná ex­kurze kolem lon­gest common sub­sequence (LCS)diffů, mě na­vedla k úvahám o rych­lej­ším LCS přes apro­xi­maci a odhady, pak jsem se dostal k rych­lým dy­na­mic time war­ping al­go­ritmům, které by šly adap­to­vat pro LCS, jen kdyby bylo možné zprů­mě­ro­vat/od­had­nout délku LCS. Na­padla mě trans­for­mace do mno­žiny před­chá­ze­jí­cích dvojic, na které by se dala apli­ko­vat Jac­car­dova po­dob­nost a tedy i minha­shinglo­ca­lity sensi­tive ha­shing (LSH). Pak jsem začal pátrat po přímém LSH pro LCS a na­ra­zil na hro­madu paperů po­pi­su­jí­cích LSH pro li­bo­volné me­t­rické pro­story, kde vzdá­le­nost je black-box funkce, nebo do­konce i pro ne­me­t­rické po­story. To všechno spus­tila jedna myš­lenka co kdyby. Začal jsem s úva­hami o efek­tivní kom­pri­maci ur­či­tých dat a skon­čil čtením o in­de­xo­vání per­mu­ta­cemi, ná­hod­ných bi­s­ek­to­rech, NAPPu, cross-po­ly­tope LSH, multi probe LSHna­vi­ga­ble small world gra­fech. To je podle mě příliš hlu­boko.

V po­dob­ném duchu jsem šel příliš hlu­boko s ná­so­be­ním matic. Díval jsem se na před­nášku Clo­jure is Not Afraid of the GPU a když před­ná­še­jící uvedl pár číslel o výkonu, začalo mě svrbět. „To dám, tyhle trhnu,“ po­mys­lel jsem si. Už jednou se mi něco ta­ko­vého trochu po­vedlo, tak proč ne teď? Stopnul jsem video, ote­vřel vim, vy­tvo­řil soubor matrix.c a začal psát.

Cesta za ho­ří­cím kře­mí­kem byla na­ko­nec úspěšná, pro­tože jsem se ±dostal na úroveň sou­čas­ného state of the art (v podobě kni­jovny ATLAS). To ale ne­stojí za po­zor­nost, pro­tože to může udělat každý. Za­jí­mavá byla cesta, kterou se člověk musí vydat, aby z pro­ce­soru dostal úplné ma­xi­mum.

Za­stávky byly ná­sle­du­jící:

první den:
naivní C (@3.2GHz)      12.000 sec
vektorizované C          1.660 sec
+ software pipelining    0.608 sec
+ blokování              0.556 sec
+ open heart surgery     0.437 sec
+ heirarchické bloky     0.310 sec
+ openmp 4 vlánka        0.080 sec

druhý a třetí den:
+ vyladění parametrů     0.290 sec
+ rozhození řádků        0.275 sec
+ nesouměrné bloky       0.260 sec
+ openmp 4 vlákna        0.067 sec

ATLAS (@3.6GHz)          0.288 sec

Jak je vidět z čísel, v jednom vlákně se dá do­sáh­nout 46× zrych­lení a pa­ra­le­li­zace je téměř per­fektní.


1) na­ï­veté

Prví verze, kterou jsem po­slepu vy­střihl za dvacet vteřin, bylo naivní ná­so­bení řádků/sloupců matic. Pro zjed­no­du­šení uva­žuji, že obě matice jsou čtver­cové a druhá matice je trans­po­no­vaná.

float dot(float *a, float *b, int len) {
  float res = 0.0f;
  for (int i = 0; i < len; i++) {
    res += a[i] * b[i];
  }
  return res;
}

Tuhle funkci je třeba spo­čí­tat pro každý prvek vý­sledné matice. Funkci je třeba za­vo­lat n2 krát a pro­vede O(n) práce. Vý­sle­dek tedy běží v O(n3) čase a pro matice o ve­li­kosti 2048x2048 to dává 12 vteřin na mém Haswellu.

for (int i = 0; i < len; i++) {
  for (int j = 0; j < len; j++) {
    res[i*len+j] = dot(a+(i*len), b+(j*len), len);
  }
}

To je dost dlouho, ale na­štěstí to může být lepší, mnohem lepší.

2) vek­to­ri­zo­vané ná­so­bení

Když chci ak­ce­le­ro­vat kód, který dělá number crun­ching, první krok musí vést k SIMD in­struk­cím. SIMD (neboli single in­struction mul­tiple data) ope­rují na ně­ko­lika po­lož­kách na­jed­nou. Ty­pic­kým pří­kla­dem je SIMD vek­to­ri­zo­vané ná­so­bení (malá ter­mi­no­lo­gická po­známka: pro řádky/sloupce matic budu po­u­ží­vat slovo „vektor“ a pro vek­tory ulo­žené v SIMD re­gis­trech pro­ce­soru budu po­u­ží­vat termín „SIMD vektor“), které vypadá ná­sle­dovně:

v₁: a₀ a₁ a₂ a₃
    *  *  *  *
v₂: b₀ b₁ b₂ b₃

v₁ a v₂ jsou zdro­jové SIMD vek­tory ulo­žené ve vek­to­ro­vých re­gis­trech a in­strukce vmulsd vy­ná­sobí čtyři od­po­ví­da­jící po­ložky pa­ra­lelně.8

Intelí CPU od dob Sandy Bridge pod­po­rují AVX/AVX2 SIMD, které má šířku 256 bitů. Do jed­noho vek­to­ro­vého re­gis­tru se vejde 8 floatů a jedna vmulsd in­strukce pa­ra­lelně pro­vede 8 ná­so­bení. Ope­race má la­tenci 5 taktů a můj Haswell dokáže každý takt od­pá­lit dvě na­jed­nou. Navíc s roz­ší­ře­ním FMA (fused mul­tiply-add) může v jedné SIMD in­strukci pro­vést vek­to­ro­vou ope­raci od­po­ví­da­jící a = b * c + d, což je přesně to, co v tomto pří­padě po­tře­buji.

Kód pro SIMD vek­to­ri­zo­vané ná­so­bení ex­pli­citně po­u­ží­va­jící in­trin­sic funkce vypadá takhle:

float dot_vec(float *a, float *b, int len) {
  // akumulátor
  __m256 d = _mm256_set1_ps(0.0);

  // vektorové součty
  for (int i = 0; i < len; i+=8) {
    __m256 aa = _mm256_load_ps(a + i);
    __m256 bb = _mm256_load_ps(b + i);
    d = _mm256_add_ps(d, _mm256_mul_ps(aa, bb));
  }

  // horizontální součet SIMD vektoru
  d = _mm256_hadd_ps(d, d);
  d = _mm256_hadd_ps(d, d);
  return ((float*)&d)[0] + ((float*)&d)[4];
}

V kódu ne­po­u­ží­vám přímo FMA in­trin­sics, ale dílčí ope­race ná­so­bení a sčí­tání. GCC je pře­vede na FMA a pře­kva­pivě vy­ge­ne­ruje lepší kód než, kdy­bych použil přímo FMA. #inulol1 6

Na konci je po­třeba udělat ho­ri­zon­tální součet po­lo­žek v SIMD vek­toru, abych ho zre­du­ko­val na jeden skalár. I pro tohle x86 po­sky­tuje SIMD in­strukci vhad­dps.

3) soft­ware pi­pe­li­ning

Další změnu, kterou jsem pro­vedl, bylo soft­wa­rové pi­pe­li­no­vání. Místo, abych ná­so­bil jeden řád­kový/sloup­cový vektor za druhým, napsal jsem kód, který na­jed­nou vy­ná­sobí dva se dvěma.

struct float4 dot_cross_vec(float *a, float *b, float *c, float *d, int len) {

  __m256 ac = _mm256_set1_ps(0.0);
  __m256 ad = _mm256_set1_ps(0.0);
  __m256 bc = _mm256_set1_ps(0.0);
  __m256 bd = _mm256_set1_ps(0.0);

  for (int i = 0; i < len; i+=8) {
    __m256 aa = _mm256_load_ps(a + i);
    __m256 bb = _mm256_load_ps(b + i);
    __m256 cc = _mm256_load_ps(c + i);
    __m256 bd = _mm256_load_ps(d + i);

    ac = _mm256_add_ps(ac, _mm256_mul_ps(aa, cc));
    ad = _mm256_add_ps(ad, _mm256_mul_ps(aa, bd));
    bc = _mm256_add_ps(bc, _mm256_mul_ps(bb, cc));
    bd = _mm256_add_ps(bd, _mm256_mul_ps(bb, bd));
  }

  struct float4 res = { hadd(ac), hadd(ad), hadd(bc), hadd(bd) };
  return res;
}

Myš­lenka je taková, že tohle povede k větším blokům kódu, které ob­sa­hují více ne­zá­vis­lých ope­rací v každé ite­raci, a hard­ware tyto ne­zá­vislé ope­race může pro­vést pa­ra­lelně. Jako bonus je po­třeba načíst jen 4 SIMD vek­tory pro 4 ná­so­bení, na­místo dvou pro jedno ná­so­bení, jak je tomu ve funkci dot_vec.2

SW pi­pe­li­no­vání sku­tečně vede ke sní­žení počtu pro­ve­de­ných in­strukcí, ale to není jeho hlavní efekt. Jde o formu bloc­kingu/ti­lingu – z paměti načtu blok dat, který pak opa­ko­vaně po­u­ží­vám. Tady je to pro blok 2x2 a zno­vupo­u­žití pro­bíhá v SIMD re­gis­trech o délce 8 floatů. Blo­ko­vání ob­vykle pro­bí­háh pro větší bloky a zno­vupo­u­žití je cílené na CPU cache. Tato změna efek­tivně se­škrtne pa­mě­ťové pře­nosy na po­lo­vinu a zá­ro­veň zvýší lo­ka­litu.

Z čísel je vidět, že SW pi­pe­li­no­vání fun­guje, pro­tože to vede k troj­ná­sob­nému zrych­lení.

4) bloc­king

Když jsem kód pro­fi­lo­val perfem, vy­ka­zo­val velké množ­ství L1-icache-load-missesLLC-load-misses udá­lostí. To zna­mená jediné – pro­gram je možná vý­po­četně efek­tivní, ale lo­ka­lita re­fe­rence není nijak slavná. Pro­gram by rád něco po­čí­tal, ale nemá co, pro­tože musí čekat na data než dorazí z paměti nebo L3 cache.

Pro­blém je v tom, že když ná­so­bím řád­kové vek­tory ai s bj, tak aspoň jeden z oněch dvou vek­torů jsem zatím ne­vi­děl, není v cache a proto ho musím stre­a­mo­vat z paměti. Sek­venční prů­chod je re­la­tivně efek­tivní, pro­tože ho CPU de­te­kuje a začne data na­čí­tat do­předu. Není ale zda­leka tak efek­tivní, jako pře­nosy přímo z L1 cache, pro­tože je li­mi­to­ván la­tencí a pro­pust­ností RAM a na za­čátku ite­race může utrpět drahý cache miss.3

Navíc pokud je ma­ti­cový vektor dlouhý a ne­ve­jde se do L1 cache, musím ho vždycky lovit z nižší úrovně cache, pro­tože prvky na za­čátku jsou v době, kdy je opět po­tře­buji, zcela jistě vy­ho­zeny z cache. Např. 8k prvků od­po­vídá 32kB paměti, což je přesně ve­li­kost L1 cache v mém stroji. Při ná­so­bení ale po­tře­buji dva takové vek­tory. O L1 cache tedy sou­peří 64kB dat s tím, že nej­starší po­ložky jsou z cache vy­ho­zeny jako první4 . V tomto pří­padě L1/L2 cache vůbec ne­vy­u­ží­vám a všechno stre­a­muji z L2/L3/RAM a to má k ide­ální si­tu­aci velice daleko.

V této si­tu­aci je­di­ným lo­gic­kým krokem je pl­no­hod­notné blo­ko­vání.

memset(res, 0, len*len*sizeof(float));

for (int tilei = 0; tilei < len; tilei += tilesize) {
  for (int tilej = 0; tilej < len; tilej += tilesize) {

    for (int block = 0; block < len; block += blocksize) {

      for (int i = tilei; i < tilei+tilesize; i++) {
        for (int j = tilej; j < tilej+tilesize; j++) {

          // dot_vec nebo sw pipelinovaná verze
          res[i*len+j] += dot_vec(a+i*len+block, b+j*len+block, blocksize)

        }
      }

    }

  }
}

Pa­ra­me­try tilesizeblocksize by měly být zvo­leny tak, aby se 2 * tilesize * blocksize * sizeof(float) vešlo do nějaké úrovně cache.

Použil jsem tilesize = 32blocksize = 128. To do­hro­mady dává 32kB da­ta­set, který se jen tak tak vejde do L1 a po­ho­dlně bude sedět v L2 (která má 256kB). Vý­sle­dek běží o něco málo rych­leji, ale ne nijak zá­sadně. Důvod je za­jí­mavý a vy­svět­lím ho o pár od­stavců níž.

5) open heart sur­gery

Blo­ko­vání vedlo ke zrych­lení, ale zá­ro­veň si vy­nu­tilo opa­ko­vané ho­ri­zon­tální součty. Ty jsou, jako všechno, celkem rychlé, ale nejsou za­darmo.

Ho­ri­zon­tální součty je možné eli­mi­no­vat za cenu oš­k­li­vého kódu. Pro každý blok/tile si budu pa­ma­to­vat me­zi­součty ne jako float* ale jako __m256*, tedy nikoli jako pole floatů, ale pole SIMD re­gis­trů. Tak můžu aku­mu­lo­vat vek­to­rové součty bez ho­ri­zon­tální re­dukce a na­ko­nec udělat jen jeden ho­ri­zon­tální součet.

Zjed­no­du­šená verze bez SW pi­pe­line:

memset(res, 0, len*len*sizeof(float));

for (int tilei = 0; tilei < len; tilei += tilesize) {
  for (int tilej = 0; tilej < len; tilej += tilesize) {

    __m256 sums[tilesize*tilesize];
    memset(sums, 0, tilesize*tilesize*sizeof(__m256));

    for (int block = 0; block < len; block += blocksize) {

      for (int i = tilei; i < tilei+tilesize; i++) {
        for (int j = tilej; j < tilej+tilesize; j++) {

          for (int i = 0; i < len; i+=8) {
            __m256 aa = _mm256_load_ps(a+i*len+block + i);
            __m256 bb = _mm256_load_ps(b + i);
            sums[i*tilesize+j] = _mm256_add_ps(d, _mm256_mul_ps(aa, bb));
          }

        }
      }

    }

    // závěrečné horizontální součty
    for (int i = tilei; i < tilei+tilesize; i++) {
      for (int j = tilej; j < tilej+tilesize; j++) {
        res[i*len+j] = hadd(i*tilesize+j)
      }
    }

  }
}

Eli­mi­nace ho­ri­zon­tál­ních součtů je pře­kva­pivě efek­tivní a vede ke zrych­lení o 20%.

6) hi­e­rar­chické bloky

Dalším lo­gic­kým krokem je hi­e­rar­chické blo­ko­vání. Místo bloků, které jsou uzpů­so­beny pro jednu úroveň cache, je lepší mít ně­ko­lik úrovní bloků, jednu pro každou úroveň cache. Jeden L1 blok ite­ruje pouze v L1 cache. Když do­končí práci a posune se na ve­d­lejší L1 blok, data budou v L2 cache. Stejně tak, když se posunu na ná­sle­du­jící L2 blok, po­třebná data budou v L3 cache.

Vý­sledný kód je oš­k­livý jako hřích. Má devět vno­ře­ných smyček, pokud dělám tří­úrov­ňové blo­ko­vání v osách x a y a ploché v ose z, nebo 12 smyček, když blo­kuji hi­e­rar­chicky ve všech osách.

Hi­e­rar­chické bloky ve vše osách by měly být rych­lejší, pro­tože se cache vy­u­žívá ide­álně, ale v reálu jsou o něco po­ma­lejší. Ne­vy­ka­zují skoro žádné L1 cache miss, ale mají hodně L3 cache miss a podle mě za to může pre­fetch, který na pozadí načítá nové a nové bloky a stíhá pro­tože CPU musí udělat řádově ku­bické množ­ství vý­po­čtů, ale po­tře­buje načíst jen kva­d­ra­tické množ­ství paměti.

Nejmenší množ­ství paměti, které je třeba načíst, by mělo být 2fn3 / √(s/(2f)), kde n je ve­li­kost matice, s ve­li­kost cache v baj­tech a f ve­li­kost floatu/double. Když to pře­po­čí­tám na L3 cache-miss, můj pro­gram jich vy­ka­zuje méně něž tohle te­o­re­tické mi­ni­mum. To uka­zuje na práci pre­fet­cheru na pozadí.

Pokud bych chtěl být moc nóbl (a cache ob­li­vi­ous) pro­chá­zel bych maticí ve stylu pro­stor vy­pl­ňu­jí­cích křivek a Mor­to­nova roz­kladu (Z-order).

7) pa­ra­lelní bu­douc­nost

Po­sled­ním krokem byla pa­ra­le­li­zace. Pře­kva­pivě to byla nej­snazší úprava, mnohem jed­no­dušší než každý z před­cho­zích kroků. Byla do­konce tak tri­vi­ální, že ne­po­tře­bo­vala žádnou změnu zdro­jo­vého kódu. Sta­čilo přidat pragmu před nej­vyšší smyčku

#pragma omp parallel for

a zkom­pi­lo­vat pro­gram s pře­pí­na­čem -fopenmp a všechno jelo téměř 4× rych­leji na 4 já­drech.

Pa­ra­le­li­zace je jed­no­du­chá, pro­tože ite­race smyčky na sobě jsou ne­zá­vislé a ne­pro­bíhá mezi nimi žádná ko­mu­ni­kace. Bez ko­mu­ni­kace, ať už ex­pli­citní pro­střed­nic­tvím zpráv nebo im­pli­citní pro­střed­nic­tvím sdí­lené paměti, je pa­ra­le­lis­mus tri­vi­ální.5 Kon­cepčně každá ite­race vy­ná­sobí jeden řád­kový vektor s jedním sloup­co­vým a vý­sle­dek zapíše do při­pra­ve­ného pole. Jediná ko­mu­ni­kace pro­bíhá na úrovni ne­vi­di­telné pro­gra­má­to­rovi – na­pří­klad false sha­ring jed­not­li­vých cache line ke kterým při­stu­puje více vláken. Ale to je jen te­o­re­tická even­tu­a­lita, pro­tože blo­ko­vání za­rov­nané na ná­so­bek cache line to zcela od­straní.

Ještě bych měl zmínit jednu věc:

Pa­ra­le­li­zace po­sky­tuje téměř per­fektní zrych­lení proto, že všechna data čte pře­vážně z cache a ne­po­tře­buje skoro nic stre­a­mo­vat z RAM. Hlavní paměť je sdí­lený zdroj s ome­ze­nou pro­pust­ností a kdyby se o něj prala 4 vlákna, mohl by se vy­čer­pat, což by vedlo k úzkému hrdlu pa­ra­le­li­zace.


I když vý­sledky vy­pa­dají pěkně, je třeba mít na paměti, že GPU jsou rych­lejší. V před­nášce, která to všechno začalo, se říká, že jeden kon­krétní model gra­fické karty je 60× rych­lejší než kon­krétní CPU. Gra­fické karty jsou pa­ra­lelní vý­po­četní mon­stra, která mají na palubě my­ri­ády pro­ce­sorů a jader a velice rych­lou paměť. Každá vý­po­četní jed­notka je sama o sobě pomalá a nehodí se pro obecné pro­gra­mo­vání, ale pro GPU platí, že síla je ve vel­kých po­čtech a pra­vi­del­nosti.

Ma­ti­cové ná­so­bení v jednom vlákně dává kolem 55 Gflops. Te­o­re­tické ma­xi­mum kaž­dého jádra mého pro­ce­soru je přitom 102 Gflops. CPU běží na 3.2 GHz, každý takt může od­pá­lit dvě FMA in­strukce, každá z nich dává 2 flops (ná­so­bení a sčí­tání) na 8 flo­a­tech. Takže stále tu je určité místo pro zlep­šení.

Bu­doucí pro­ce­sory od Intelu (ser­ve­rové Sky­lake a Can­non­lake) mají při­nést pod­poru pro AVX-512. To je dva­krát širší než AVX/AVX2 a tedy by mohlo vést k dvoj­ná­sob­nému zrych­lení.

Ještě po­známka: Během celé téhle pro­gra­mo­vací es­ka­pády jsem se snažil zlep­šit vy­u­žití cache ex­pli­cit­ními prefetch in­struk­cemi, ale nikdy to nikam ne­vedlo. To co dělá sám od sebe bylo vždy víc než do­sta­tečné.

8) po uzá­věrce

Den poté, co jsem dopsal první verzi to­ho­hle článku jsem zasedl ke stro­jům a udělal pár op­ti­ma­li­zací, které vedly k dalším drob­ným zrych­le­ním.

+ vyladění parametrů     0.290 sec (à la ATLAS)
+ rozhození řádků        0.275 sec
+ nesouměrné bloky       0.260 sec
+ openmp 4 vlákna        0.067 sec

Prvním z nich bylo au­to­ma­tické vy­la­dění pa­ra­me­trů – těch je tu celá řada: např. ve­li­kosti tří úrovní bloků (tile_size), kolik po­lo­žek řádku matice se načte pro každý blok (block_size) nebo různé pře­pí­nače a volby pro GCC. Ně­které z nich vedou ke zrych­lení, jiné k pro­padu výkonu. Nemá smysl je všechny zkou­šet ma­nu­álně, ale nej­lepší kom­bi­naci hledat au­to­ma­ticky. Uká­zalo se, že nej­rych­lejší je dvou­úrov­ňová struk­tura: bloky ve­li­kosti 8 pro L1, 256 pro L2/L3 a 32 floatů na­jed­nou.

To mi bylo divné, pro­tože L1 blok zabere jen 8324B = 1kB paměti. Čekal bych, že lep­šího výkonu se do­sáhne, když L1 blok vy­u­žije víc z 32kB L1 cache. Ale na­vý­šení L1 tile_size na 16 vedlo ke zpo­ma­lení o 30% a ná­růstu počtu L1 cache-miss z 13.25% na 34%.

Cache je or­ga­ni­zo­vaná takto:

Jde o 32kB hash-ta­bulku im­ple­men­to­va­nou v hard­ware tvo­ře­nou 64B bloky cache-line. Každý bucket (set) ta­bulky má fixní ve­li­kost/aso­ci­a­ti­vitu (v tomto pří­padě 8, tedy 8-way). Ha­sho­vání je za­lo­žené na adrese daného slova/bajtu/če­ho­koli: (adresa >> 6) & ((1 << 6) - 1). Vy­sekne pár níz­kých bitů. Určité slovo tedy může být jen v jedné mno­žině, která má 8 pozic. Když jsou všechny pozice za­pl­něné, musím ně­kte­rou vy­kli­dit.

A to je právě ono. Když jsou dvě po­ložky od sebe vzdá­lené o ná­so­bek 4kB, spad­nou do stej­ného bucketu v L1. Moje tes­to­vací matice má ve­li­kost 2048x2084 a řádky, které mají délku 8kB, jsou na­sklá­dané v paměti jeden za druhým. Proto k-tý prvek z jed­noho řádku je o 8kB od­sa­zen od k-tého prvku z ná­sle­du­jí­cího řádku. Když tedy zvolím L1 tile_size = 16, pra­cuji s 16 slovy, které spad­nou do stej­ného bucketu v cache, ale ten má aso­ci­a­ti­vitu pouze 8. Tím pádem se ak­tivní da­ta­set ne­u­stále vy­ha­zuje z L1.

Řešení je jed­no­du­ché – stačí řádky trochu od­sa­dit, aby nebyly za­rov­nané na moc­ninu dvou. Když každý řádek vycpu 64 bajty, cache line s prv­ními prvky ne­spad­nou do stej­ného bucketu, ale roz­pro­střou se po celé L1 cache.7

Tato změna sní­žila počet L1 cache miss z 13.25% na 12.5% a ubrala ~15 ms z běhu pro­gramu.

Teď bych doufal, že na­vý­šení L1 bloku z 8 na 16 ko­nečně povede ke zrych­lení. Ale není tomu tak. S L1 tile_size=16 pro­gram běží po­ma­leji a perf uka­zuje zna­telný nárůst počtu vy­ko­na­ných in­strukcí. To může zna­me­nat že vnitřní smyčka je na­jed­nou příliš dlouhá a GCC se ji nej­spíš roz­hodl ne­odro­lo­vat a zkom­pi­lo­val ji jinak, po­ma­leji.

Jestliže tomu tak je, stačí zavést ne­sou­měrné bloky – vnitřní smyčka má délku 8, aby GCC byl spo­ko­jený, a vnější má délku 16, abych lépe využil L1 cache. Sku­tečně to fun­go­valo, podíl L1 cache miss spadl z 12.5% na 10.25%, pro­gram běží o dal­ších ~15 ms rych­leji a dává kolem 66 Gflops v jednom vlákně.

A co víc: pa­ra­le­li­zace je pořád dobrá. Na čtyřech vlák­nech běží 3.88× rych­leji.

O blo­ko­vání a aso­ci­a­ti­vitě cache se píše tady.

9) jen pro ty, kteří se nebojí ris­ko­vat

To že jsem dosáhl parity s kni­hov­nou, kterou na­psali lidé, kteří nejsou di­le­tanti jako já, je sice pěkné, ale co dál? Je možné ještě zrych­lit bez změny hard­ware?

Jedna z mož­ností, jak zrych­lit vý­po­čet a vy­ma­nit se o okovů O(n3) je začít od­ha­do­vat. Když oželím přes­nost, můžu vý­sle­dek apro­xi­mo­vat. Zběžně jsem na­hlédl do li­te­ra­tury, abych zjis­til co a jak. Mož­ností je ne­pře­berně, od jed­no­du­chého vzor­ko­vání až po ra­fi­no­va­nosti, které mi jdou přes hlavu.

Jedna z nich je třeba použít něco, co už znám: ske­čo­vání a LSH. Pomocí ná­hod­ných nadro­vin nebo cross-po­ly­tope můžu získat odhad pro ko­si­no­vou po­dob­nost

cossim(a, b) = a . b / (|a| * |b|)

Když tedy nor­ma­li­zuji řádky jedné matice a sloupce druhé, vy­po­čí­tám z nich skeče kratší než pů­vodní vek­tory a tyto skeče pak po­u­žiji k odhadu sou­činu nor­ma­li­zo­va­ných vektor a vý­sled­nou hod­notu vy­ná­so­bím dél­kami pů­vod­ních vek­torů do­stanu odhad ma­ti­co­vého ná­so­bení. Odhad je velice rychlý, pro­tože počítá Ha­m­min­govu vzdá­le­nost na bi­to­vých vek­to­rech, která se dá zjis­tit velice rychle kom­bi­nací xorpopcnt in­strukcí (nebo rych­leji nebo ještě rych­leji).

Nic ta­ko­vého jsem ale ne­pro­gra­mo­val, pro­tože, jak je z délky článku patrné, už tak jsem na tomhle pro­blému strá­vil víc času, než jsem chtěl.


Pod čarou:


Dále k tématu:


Pozn:

  1. Po­dob­ného vý­sledku by se možná dalo do­sáh­nout au­to­ma­ticky, kdy­bych se snažil pře­svěd­čit GCC, aby smyčku vek­to­ri­zo­val.
  2. Tech­nicky vzato by GCC mohl sám dojít k něčemu po­dob­nému, kdyby byl do­sta­tečně od­vážný, trans­po­no­val vno­řené smyčky, odro­lo­val je, věřil by si se za­rov­ná­ním paměti a počtem ite­rací a tak po­dobně. Pak by dostal velký blok vmovapsvfmaddXXXps in­strukcí, které by mohl velice efek­tivně na­plá­no­vat a eli­mi­no­vat du­pli­citní vmovaps. Tech­nicky by to asi šlo, ale sám nevím jak ho k tomu pře­svěd­čit.
  3. Pro­pust­nost RAM je za­ne­dba­telná v malá v po­rov­nání s L1/L2 cache. Nej­rych­lejší čtyř­ka­ná­lová DDR4 dokáže pro­táh­nout ±60GB/s s la­tencí 60 ns, DDR3 v mém stroji dají něco kolem 20GB/s, L1 každý tak pro­tlačí 2x32N, což dává do­hro­mady 200 GB/s na jedno jádro s la­tencí 3 takty a 800GB/s na všechna čtyři jádra.
  4. Cache In­te­lích pro­ce­sorů pro eviction policy po­u­žívá apro­xi­maci LRU pro L1/L2 a adap­tivní tech­niku pro L3, která se snaží vy­po­řá­dat právě s pří­pa­dem, kdy pro­jedu data a na každý prvek sáhnu jen jednou.
  5. viz jeden slajd v Grim hard­ware re­a­li­ties of functi­o­nal pro­gra­m­ming
  6. Na ně­kte­rých stro­jích pře­chod z tříad­re­so­vých FMA do dvou­ad­re­so­vých AVX in­strukcí obnáší malou pe­na­li­zaci. Proto může být rych­lejší použít čistě jen FMA nebo jen AVX in­strukce, nebo vy­ko­nat ně­ko­lik AVX in­strukcí a pak ně­ko­lik FMA in­strukcí a pe­na­li­zaci tak amor­ti­zo­vat.
  7. Více k tomuto tématu je v paperu Array la­y­outs for com­pa­ri­son-based sear­chingdal­ších člán­cích.
  8. SIMD in­strukce můžou najít ne­če­kané vy­u­žití, jako na­pří­klad: Fast Sorted-Set In­ter­section using SIMD In­structi­ons.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články