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

Mám tendenci jít vždycky příliš hluboko.

Například moje nedávná exkurze kolem longest common subsequence (LCS) a diffů, mě navedla k úvahám o rychlejším LCS přes aproximaci a odhady, pak jsem se dostal k rychlým dynamic time warping algoritmům, které by šly adaptovat pro LCS, jen kdyby bylo možné zprůměrovat/odhadnout délku LCS. Napadla mě transformace do množiny předcházejících dvojic, na které by se dala aplikovat Jaccardova podobnost a tedy i minhashing a locality sensitive hashing (LSH). Pak jsem začal pátrat po přímém LSH pro LCS a narazil na hromadu paperů popisujících LSH pro libovolné metrické prostory, kde vzdálenost je black-box funkce, nebo dokonce i pro nemetrické postory. To všechno spustila jedna myšlenka co kdyby. Začal jsem s úvahami o efektivní komprimaci určitých dat a skončil čtením o indexování permutacemi, náhodných bisektorech, NAPPu, cross-polytope LSH, multi probe LSH a navigable small world grafech. To je podle mě příliš hluboko.

V podobném duchu jsem šel příliš hluboko s násobením matic. Díval jsem se na přednášku Clojure is Not Afraid of the GPU a když přednášející uvedl pár číslel o výkonu, začalo mě svrbět. "To dám, tyhle trhnu," pomyslel jsem si. Už jednou se mi něco takového trochu povedlo, tak proč ne teď? Stopnul jsem video, otevřel vim, vytvořil soubor matrix.c a začal psát.

Cesta za hořícím křemíkem byla nakonec úspěšná, protože jsem se ±dostal na úroveň současného state of the art (v podobě knijovny ATLAS). To ale nestojí za pozornost, protože to může udělat každý. Zajímavá byla cesta, kterou se člověk musí vydat, aby z procesoru dostal úplné maximum.

Zastávky byly následují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á dosáhnout 46x zrychlení a paralelizace je téměř perfektní.


1) naïveté

Prví verze, kterou jsem poslepu vystřihl za dvacet vteřin, bylo naivní násobení řádků/sloupců matic. Pro zjednodušení uvažuji, že obě matice jsou čtvercové a druhá matice je transponovaná.

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 zavolat n2 krát a provede O(n) práce. Výsledek tedy běží v O(n3) čase a pro matice o velikosti 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) vektorizované násobení

Když chci akcelerovat kód, který dělá number crunching, první krok musí vést k SIMD instrukcím. SIMD (neboli single instruction multiple data) operují na několika položkách najednou. Typickým příkladem je SIMD vektorizované násobení (malá terminologická poznámka: pro řádky/sloupce matic budu používat slovo "vektor" a pro vektory uložené v SIMD registrech procesoru budu používat termín "SIMD vektor"), které vypadá následovně:

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

v₁ a v₂ jsou zdrojové SIMD vektory uložené ve vektorových registrech a instrukce vmulsd vynásobí čtyři odpovídající položky paralelně.8

Intelí CPU od dob Sandy Bridge podporují AVX/AVX2 SIMD, které má šířku 256 bitů. Do jednoho vektorového registru se vejde 8 floatů a jedna vmulsd instrukce paralelně provede 8 násobení. Operace má latenci 5 taktů a můj Haswell dokáže každý takt odpálit dvě najednou. Navíc s rozšířením FMA (fused multiply-add) může v jedné SIMD instrukci provést vektorovou operaci odpovídající a = b * c + d, což je přesně to, co v tomto případě potřebuji.

Kód pro SIMD vektorizované násobení explicitně používající intrinsic 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 nepoužívám přímo FMA intrinsics, ale dílčí operace násobení a sčítání. GCC je převede na FMA a překvapivě vygeneruje lepší kód než, kdybych použil přímo FMA. #inulol1 6

Na konci je potřeba udělat horizontální součet položek v SIMD vektoru, abych ho zredukoval na jeden skalár. I pro tohle x86 poskytuje SIMD instrukci vhaddps.

3) software pipelining

Další změnu, kterou jsem provedl, bylo softwarové pipelinování. Místo, abych násobil jeden řádkový/sloupcový vektor za druhým, napsal jsem kód, který najednou vyná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é obsahují více nezávislých operací v každé iteraci, a hardware tyto nezávislé operace může provést paralelně. Jako bonus je potřeba načíst jen 4 SIMD vektory pro 4 násobení, namísto dvou pro jedno násobení, jak je tomu ve funkci dot_vec.2

SW pipelinování skutečně vede ke snížení počtu provedených instrukcí, ale to není jeho hlavní efekt. Jde o formu blockingu/tilingu - z paměti načtu blok dat, který pak opakovaně používám. Tady je to pro blok 2x2 a znovupoužití probíhá v SIMD registrech o délce 8 floatů. Blokování obvykle probíháh pro větší bloky a znovupoužití je cílené na CPU cache. Tato změna efektivně seškrtne paměťové přenosy na polovinu a zároveň zvýší lokalitu.

Z čísel je vidět, že SW pipelinování funguje, protože to vede k trojnásobnému zrychlení.

4) blocking

Když jsem kód profiloval perfem, vykazoval velké množství L1-icache-load-misses a LLC-load-misses událostí. To znamená jediné - program je možná výpočetně efektivní, ale lokalita reference není nijak slavná. Program by rád něco počítal, ale nemá co, protože musí čekat na data než dorazí z paměti nebo L3 cache.

Problém je v tom, že když násobím řádkové vektory ai s bj, tak aspoň jeden z oněch dvou vektorů jsem zatím neviděl, není v cache a proto ho musím streamovat z paměti. Sekvenční průchod je relativně efektivní, protože ho CPU detekuje a začne data načítat dopředu. Není ale zdaleka tak efektivní, jako přenosy přímo z L1 cache, protože je limitován latencí a propustností RAM a na začátku iterace může utrpět drahý cache miss.3

Navíc pokud je maticový vektor dlouhý a nevejde se do L1 cache, musím ho vždycky lovit z nižší úrovně cache, protože prvky na začátku jsou v době, kdy je opět potřebuji, zcela jistě vyhozeny z cache. Např. 8k prvků odpovídá 32kB paměti, což je přesně velikost L1 cache v mém stroji. Při násobení ale potřebuji dva takové vektory. O L1 cache tedy soupeří 64kB dat s tím, že nejstarší položky jsou z cache vyhozeny jako první4 . V tomto případě L1/L2 cache vůbec nevyužívám a všechno streamuji z L2/L3/RAM a to má k ideální situaci velice daleko.

V této situaci jediným logickým krokem je plnohodnotné bloková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)

        }
      }

    }

  }
}

Parametry tilesize a blocksize by měly být zvoleny tak, aby se 2 * tilesize * blocksize * sizeof(float) vešlo do nějaké úrovně cache.

Použil jsem tilesize = 32 a blocksize = 128. To dohromady dává 32kB dataset, který se jen tak tak vejde do L1 a pohodlně bude sedět v L2 (která má 256kB). Výsledek běží o něco málo rychleji, ale ne nijak zásadně. Důvod je zajímavý a vysvětlím ho o pár odstavců níž.

5) open heart surgery

Blokování vedlo ke zrychlení, ale zároveň si vynutilo opakované horizontální součty. Ty jsou, jako všechno, celkem rychlé, ale nejsou zadarmo.

Horizontální součty je možné eliminovat za cenu ošklivého kódu. Pro každý blok/tile si budu pamatovat mezisoučty ne jako float* ale jako __m256*, tedy nikoli jako pole floatů, ale pole SIMD registrů. Tak můžu akumulovat vektorové součty bez horizontální redukce a nakonec udělat jen jeden horizontální součet.

Zjednodušená verze bez SW pipeline:

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

  }
}

Eliminace horizontálních součtů je překvapivě efektivní a vede ke zrychlení o 20%.

6) hierarchické bloky

Dalším logickým krokem je hierarchické blokování. Místo bloků, které jsou uzpůsobeny pro jednu úroveň cache, je lepší mít několik úrovní bloků, jednu pro každou úroveň cache. Jeden L1 blok iteruje pouze v L1 cache. Když dokončí práci a posune se na vedlejší L1 blok, data budou v L2 cache. Stejně tak, když se posunu na následující L2 blok, potřebná data budou v L3 cache.

Výsledný kód je ošklivý jako hřích. Má devět vnořených smyček, pokud dělám tříúrovňové blokování v osách x a y a ploché v ose z, nebo 12 smyček, když blokuji hierarchicky ve všech osách.

Hierarchické bloky ve vše osách by měly být rychlejší, protože se cache využívá ideálně, ale v reálu jsou o něco pomalejší. Nevykazují skoro žádné L1 cache miss, ale mají hodně L3 cache miss a podle mě za to může prefetch, který na pozadí načítá nové a nové bloky a stíhá protože CPU musí udělat řádově kubické množství výpočtů, ale potřebuje načíst jen kvadratické 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 velikost matice, s velikost cache v bajtech a f velikost floatu/double. Když to přepočítám na L3 cache-miss, můj program jich vykazuje méně něž tohle teoretické minimum. To ukazuje na práci prefetcheru na pozadí.

Pokud bych chtěl být moc nóbl (a cache oblivious) procházel bych maticí ve stylu prostor vyplňujících křivek a Mortonova rozkladu (Z-order).

7) paralelní budoucnost

Posledním krokem byla paralelizace. Překvapivě to byla nejsnazší úprava, mnohem jednodušší než každý z předchozích kroků. Byla dokonce tak triviální, že nepotřebovala žádnou změnu zdrojového kódu. Stačilo přidat pragmu před nejvyšší smyčku

#pragma omp parallel for

a zkompilovat program s přepínačem -fopenmp a všechno jelo téměř 4x rychleji na 4 jádrech.

Paralelizace je jednoduchá, protože iterace smyčky na sobě jsou nezávislé a neprobíhá mezi nimi žádná komunikace. Bez komunikace, ať už explicitní prostřednictvím zpráv nebo implicitní prostřednictvím sdílené paměti, je paralelismus triviální.5 Koncepčně každá iterace vynásobí jeden řádkový vektor s jedním sloupcovým a výsledek zapíše do připraveného pole. Jediná komunikace probíhá na úrovni neviditelné programátorovi - například false sharing jednotlivých cache line ke kterým přistupuje více vláken. Ale to je jen teoretická eventualita, protože blokování zarovnané na násobek cache line to zcela odstraní.

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

Paralelizace poskytuje téměř perfektní zrychlení proto, že všechna data čte převážně z cache a nepotřebuje skoro nic streamovat z RAM. Hlavní paměť je sdílený zdroj s omezenou propustností a kdyby se o něj prala 4 vlákna, mohl by se vyčerpat, což by vedlo k úzkému hrdlu paralelizace.


I když výsledky vypadají pěkně, je třeba mít na paměti, že GPU jsou rychlejší. V přednášce, která to všechno začalo, se říká, že jeden konkrétní model grafické karty je 60x rychlejší než konkrétní CPU. Grafické karty jsou paralelní výpočetní monstra, která mají na palubě myriády procesorů a jader a velice rychlou paměť. Každá výpočetní jednotka je sama o sobě pomalá a nehodí se pro obecné programování, ale pro GPU platí, že síla je ve velkých počtech a pravidelnosti.

Maticové násobení v jednom vlákně dává kolem 55 Gflops. Teoretické maximum každého jádra mého procesoru je přitom 102 Gflops. CPU běží na 3.2 GHz, každý takt může odpálit dvě FMA instrukce, každá z nich dává 2 flops (násobení a sčítání) na 8 floatech. Takže stále tu je určité místo pro zlepšení.

Budoucí procesory od Intelu (serverové Skylake a Cannonlake) mají přinést podporu pro AVX-512. To je dvakrát širší než AVX/AVX2 a tedy by mohlo vést k dvojnásobnému zrychlení.

Ještě poznámka: Během celé téhle programovací eskapády jsem se snažil zlepšit využití cache explicitními prefetch instrukcemi, ale nikdy to nikam nevedlo. To co dělá sám od sebe bylo vždy víc než dostatečné.

8) po uzávěrce

Den poté, co jsem dopsal první verzi tohohle článku jsem zasedl ke strojům a udělal pár optimalizací, které vedly k dalším drobným zrychlení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 automatické vyladění parametrů - těch je tu celá řada: např. velikosti tří úrovní bloků (tile_size), kolik položek řádku matice se načte pro každý blok (block_size) nebo různé přepínače a volby pro GCC. Některé z nich vedou ke zrychlení, jiné k propadu výkonu. Nemá smysl je všechny zkoušet manuálně, ale nejlepší kombinaci hledat automaticky. Ukázalo se, že nejrychlejší je dvouúrovňová struktura: bloky velikosti 8 pro L1, 256 pro L2/L3 a 32 floatů najednou.

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

Cache je organizovaná takto:

Jde o 32kB hash-tabulku implementovanou v hardware tvořenou 64B bloky cache-line. Každý bucket (set) tabulky má fixní velikost/asociativitu (v tomto případě 8, tedy 8-way). Hashování je založené na adrese daného slova/bajtu/čehokoli: (adresa >> 6) & ((1 << 6) - 1). Vysekne pár nízkých bitů. Určité slovo tedy může být jen v jedné množině, která má 8 pozic. Když jsou všechny pozice zaplněné, musím některou vyklidit.

A to je právě ono. Když jsou dvě položky od sebe vzdálené o násobek 4kB, spadnou do stejného bucketu v L1. Moje testovací matice má velikost 2048x2084 a řádky, které mají délku 8kB, jsou naskládané v paměti jeden za druhým. Proto k-tý prvek z jednoho řádku je o 8kB odsazen od k-tého prvku z následujícího řádku. Když tedy zvolím L1 tile_size = 16, pracuji s 16 slovy, které spadnou do stejného bucketu v cache, ale ten má asociativitu pouze 8. Tím pádem se aktivní dataset neustále vyhazuje z L1.

Řešení je jednoduché - stačí řádky trochu odsadit, aby nebyly zarovnané na mocninu dvou. Když každý řádek vycpu 64 bajty, cache line s prvními prvky nespadnou do stejného bucketu, ale rozprostř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 programu.

Teď bych doufal, že navýšení L1 bloku z 8 na 16 konečně povede ke zrychlení. Ale není tomu tak. S L1 tile_size=16 program běží pomaleji a perf ukazuje znatelný nárůst počtu vykonaných instrukcí. To může znamenat že vnitřní smyčka je najednou příliš dlouhá a GCC se ji nejspíš rozhodl neodrolovat a zkompiloval ji jinak, pomaleji.

Jestliže tomu tak je, stačí zavést nesouměrné bloky - vnitřní smyčka má délku 8, aby GCC byl spokojený, a vnější má délku 16, abych lépe využil L1 cache. Skutečně to fungovalo, podíl L1 cache miss spadl z 12.5% na 10.25%, program běží o dalších ~15 ms rychleji a dává kolem 66 Gflops v jednom vlákně.

A co víc: paralelizace je pořád dobrá. Na čtyřech vláknech běží 3.88x rychleji.

O blokování a asociativitě cache se píše tady.

9) jen pro ty, kteří se nebojí riskovat

To že jsem dosáhl parity s knihovnou, kterou napsali lidé, kteří nejsou diletanti jako já, je sice pěkné, ale co dál? Je možné ještě zrychlit bez změny hardware?

Jedna z možností, jak zrychlit výpočet a vymanit se o okovů O(n3) je začít odhadovat. Když oželím přesnost, můžu výsledek aproximovat. Zběžně jsem nahlédl do literatury, abych zjistil co a jak. Možností je nepřeberně, od jednoduchého vzorkování až po rafinovanosti, které mi jdou přes hlavu.

Jedna z nich je třeba použít něco, co už znám: skečování a LSH. Pomocí náhodných nadrovin nebo cross-polytope můžu získat odhad pro kosinovou podobnost

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

Když tedy normalizuji řádky jedné matice a sloupce druhé, vypočítám z nich skeče kratší než původní vektory a tyto skeče pak použiji k odhadu součinu normalizovaných vektor a výslednou hodnotu vynásobím délkami původních vektorů dostanu odhad maticového násobení. Odhad je velice rychlý, protože počítá Hammingovu vzdálenost na bitových vektorech, která se dá zjistit velice rychle kombinací xor a popcnt instrukcí (nebo rychleji nebo ještě rychleji).

Nic takového jsem ale neprogramoval, protože, jak je z délky článku patrné, už tak jsem na tomhle problému strávil víc času, než jsem chtěl.


Pod čarou:


Dále k tématu:


Pozn:

  1. Podobného výsledku by se možná dalo dosáhnout automaticky, kdybych se snažil přesvědčit GCC, aby smyčku vektorizoval.
  2. Technicky vzato by GCC mohl sám dojít k něčemu podobnému, kdyby byl dostatečně odvážný, transponoval vnořené smyčky, odroloval je, věřil by si se zarovnáním paměti a počtem iterací a tak podobně. Pak by dostal velký blok vmovaps a vfmaddXXXps instrukcí, které by mohl velice efektivně naplánovat a eliminovat duplicitní vmovaps. Technicky by to asi šlo, ale sám nevím jak ho k tomu přesvědčit.
  3. Propustnost RAM je zanedbatelná v malá v porovnání s L1/L2 cache. Nejrychlejší čtyřkanálová DDR4 dokáže protáhnout ±60GB/s s latencí 60 ns, DDR3 v mém stroji dají něco kolem 20GB/s, L1 každý tak protlačí 2x32N, což dává dohromady 200 GB/s na jedno jádro s latencí 3 takty a 800GB/s na všechna čtyři jádra.
  4. Cache Intelích procesorů pro eviction policy používá aproximaci LRU pro L1/L2 a adaptivní techniku pro L3, která se snaží vypořádat právě s případem, kdy projedu data a na každý prvek sáhnu jen jednou.
  5. viz jeden slajd v Grim hardware realities of functional programming
  6. Na některých strojích přechod z tříadresových FMA do dvouadresových AVX instrukcí obnáší malou penalizaci. Proto může být rychlejší použít čistě jen FMA nebo jen AVX instrukce, nebo vykonat několik AVX instrukcí a pak několik FMA instrukcí a penalizaci tak amortizovat.
  7. Více k tomuto tématu je v paperu Array layouts for comparison-based searching a dalších článcích.
  8. SIMD instrukce můžou najít nečekané využití, jako například: Fast Sorted-Set Intersection using SIMD Instructions.
@kaja47, kaja47@k47.cz, deadbeef.k47.cz, starší články