Dobrý odhad vydá za tisíc přesných počtů

Původně jsem chtěl napsat něco obecného o pravděpodob­nostních algoritmech a užitečnosti rychlých odhadů, ale nemám na to sílu. Místo toho jsem vybral pár článků a paperů, které se týkají tématu, roztřídil je na obecné a ty o Bloom filtrech, HyperLogLogách, Count-min skečích a další. Tady to máte:

obecné

Bloom

HyperLogLog

Count-min sketch

Locality sensitive hashing

ostatní

místní

Flattr this!

Posted in Algo, Paper | Leave a comment

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

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/od­hadnout 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 46× 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 2048×2048 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/tilin­gu – z paměti načtu blok dat, který pak opakovaně používám. Tady je to pro blok 2×2 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ěř 4× 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 Kon­cepč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 60× 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é.

7) 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/asoci­ativitu (v tomto případě 8, tedy 8-way). Hashování je založené na adrese daného slova/bajtu/če­hokoli: (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 2048×2084 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 ~15ms 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 ~15ms 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.88× rychleji.

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

8) 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:

  • Funkce násobení tady.
  • Použitím zmíněných technik (SW pipelinování a hierarchické blokování) je možné zrychlit toto z 4.2 na 2.4 vteřiny, tedy o 75%.

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í 60ns, DDR3 v mém stroji dají něco kolem 20GB/s, L1 každý tak protlačí 2×32N, 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.

Flattr this!

Posted in Uncategorized | 2 Comments

Jazyk D a radix sort

V poslední době jsem se začal učit jazyk D. Chtěl jsem poznat nějaký nový jazyk, nejlépe něco kompilovaného, co mi poskytne dobrou kontrolou nad layoutem paměti. Zkoušel jsem Rust, ale nějak moc mě nechytl, vždycky mi připadalo, že učící křivka začíná strašlivým hrbolem. Chtěl jsem začít tak, že napíšu nějaký jednoduchý program, ale v případě Rustu se musím prokousat přes reference, vypůjčování a systém vlastnictví objektů/paměti, než se do něčeho můžu pustit, a ani potom nemám dojem, že znám všechno potřebné. Chtěl bych začít stylem na hulváta: Něco o syntaxi, pár datových struktur a nějaké ty idiomy, pustit vim a začít. Ale v Rustu vždycky pochybuji: má to být objekt, reference na objekt, refcountovaná reference nebo atomicky refcountovaná reference. Co z toho mám použít a jak?

Pak jsem se pustil do D a neslo se to v duchu: Znáš C/C++? Jo? Tak je to podobné, jen tohle, tohle a tohle je jinak, má to v sobě CG, kontroluje rozsahy polí a o pointery se není třeba příliš starat (segfaultnul jsem jen jednou), takhle uděláš dynamické pole, takhle asociativní, přeji příjemné kódování. Hotovo.

D má tu příjemnou (a překvapivou) vlastnost, že všechno prostě funguje. Když nevím, jak něco udělat, něco zkusím a ejhle: Ono to je správné řešení. I když jazyk neznám, působí familiárním dojmem, má bohatou standardní knihovnu a proto je s ním jednoduché začít.

Jen typový systém mě trochu děsí. Ze Scaly jsem zvyklý na něco rigorózního a systematického. D nemá nic takového, jen nekonečné instanciace šablon. To je na jednu stranu fajn, protože můžu přesně řídit, jak a pro co se mi specializuje kód, ale na druhou stranu z toho mám pocit, jako kdybych proplouval mezihvězdným prostorem: Není tam nic stabilního, o co bych se mohl opřít.


V rámci studia jsem v D napsal radix sort (zdrojáky tady). Jde překvapivě pouze o 40 řádků kódu a i když to nezvládá znaménkové typy nebo floaty, jde o templatovaný kód pro všechny neznaménkové inty1. To není nijak zlé, není se třeba zdržovat žádnou manuální specializací, která by byla potřeba v Javě/Scale, abych pořádně rozpálil křemík.

Oproti minulé verzi, kterou jsem vyseknul ve Scale, je tohle provedení v jazyku, který poskytuje kontrolu nad layoutem paměti, zaručenou alokaci na zásobníku a možností obejít kontroly rozsahů polí, rychlejší. Ale hlavně je stabilně rychlejší a můžu se na to spolehnout. Nemusím se strachovat, aby JIT všechno devirtualizoval a inlinoval.

Čas na element vstupního pole je konstantní od maličkých polí kolem 64 elementů až po gigantické s délkou 256M. Když je délka vstupního intu konstantní, jde skutečně o lineární řazení.

Následující tabulka je pro pole 64M intů (256MB kus paměti). Porovnává radix sort s řadícím algo ze standardní knihovny D (mělo by jít o introsort nebo timsort).

typ radix sort std.algorithm­.sorting.sort
ubyte 1.5 ns/el 42 ns/el
ushort 3.7 ns/el 66 ns/el
uint 9.5 ns/el 89 ns/el
ulong 28.5 ns/el 91 ns/el

s.a.s.sort má O(n log n) složitost, je pomalejší a jeho čas na element roste s velikostí vstupu. Radix sort 256G intů se stále pohybuje pod 10ns/položku, U 1G intu to spadne na 15ns na položku. To se však dá čekat. Vstupní pole má 4GB (+ 4GB originál v testovacím kódu) a dočasný pracovní prostor zabírá další 4GB (tento radix sort není in-place2). Dohromady to zabírá většinu paměti na mém stroji, který mohl občas swapovat.

Z podrobnějšího zkoumání to vypadá, že algoritmus je limitován pamětí a/nebo cache. perf hlásí velké množství L1-dcache-load-misses a také LLC-load-misses a LLC-store-misses, ale skoro žádné dTLB-load-misses. Když program čte data lineárně, ale cache selhávají, může to znamenat, že prefetch nestíhá data sypat z paměti. LLC-*-misses skoro zmizí pro menší vstupní pole intů a čas spadne na 8.5 ns/el.

ulong je oproti uint verzi 3× pomalejší. To se dá vysvětlit tak, že jednak dělá 2× více iterací (jedna iterace pro každý byte) a navíc při každé iteraci přehazuje dvojnásobek dat a přitom dělá jen velice málo práce.

Hlavní smyčka vypadá takhle:

foreach(v; arr) {
  uint c = (v >> (b * 8)) & 0xff;
  tmparr.ptr[offsets[c]] = v;
  offsets[c] += 1;
}

Jen pár instrukcí, jedna potenciální datová závislost, dvakrát čtu z paměti, jednou do ní zapisuji, dělám 4 iterace + jednu pro spočítání globálního histogramu. Pro seřazení 4B intu musím pohnout 52B dat. A paralelizace nepřinesla vůbec žádné zrychlení. Na 4 vláknech běží pořád něco přes 10ns/el. Paralelizace přinesla určité zrychlení. Na 4 jádrech běží ~2× rychleji. Poloviční zrychlení oproti ideálu se dá čekat, protože program provede 2× tolik instrukcí (rozdělení práce a spojení výsledků není zadarmo). Tahle verze přehazuje data z/do paměti kolem tempem přes 10GB/s.

D mi připadá jako dobrý jazyk pro psaní těchto nízkoúrovňových věcí. Dovolí mi dostat z hardware maximum, ale nepůsobí neohrabaně jako C nebo vyloženě nepřátelsky jako C++.


Dále k tématu:


Pozn:

  1. D má tu příjemnou vlastnost, že jeho celočíselné typy mají stejné délky jako v Javě: byte (1B), short/ushort (2B), int/uint (4B) a long/ulong (8B)
  2. Na rozdíl od PARADISu

Flattr this!

Posted in D, low level, řazení | Leave a comment

Monoid

Jasně si vzpomínám, jak jsme kdysi dávno na vysoké škole probírali pologrupy a monoidy a já si pomyslel: „K čemu je to dobré?“ Šlo o upřímnou otázku, nikoli o znevažování všeho, co nemá okamžité uplatnění. Zajímalo mě to, protože když jsem věděl, kde se daná věc používá, dodalo mi to trochu motivace a získal jsem rámec do kterého zasadit, co jsem se naučil. Tu otázku jsem nikdy nevyslovil nahlas a proto se mi nedostalo žádné odpovědi a tak jsem žil v ignoranci než mě ze školy bez pocty vyhodili a pak (aby si byli jisti) mě vyhodili ještě jednou.

Nicméně teď už to vím. Monoidy a pologrupy jsou velice užitečné algebraické struktury a abstrakce, které se uplatní všude tam, kde se dělají agregace, které byť jenom vzdáleně vypadají jako sčítání.

Monoid je množina T s jednou asociativní binární operací a neutrálním prvkem. V kódu může vypadat třeba takto:

trait Monoid[T] {
  def zero: T
  def plus(l: T, r: T): T
}

Velice podobně je definován i v knihovně Algebird, kterou budu v tomto textu používat jako příklad.

Aby byl monoid monoidem, musí splňovat dva axiomy.

  1. Binární operace musí být asociativní. Nezáleží tedy kam dám závorky, výsledek bude vždycky stejný.
a ⊕ (b ⊕ c) = (a ⊕ b) ⊕ c
  1. Nulový prvek (identita) nemá žádný vliv ať už ho přičtu zleva nebo zprava.
x ⊕ 0 = x = 0 ⊕ x

To je celá definice monoidu, nekomplikovaná a strohá. Vypadá skoro až příliš jednoduše na to, aby byla k něčemu užitečná. Když se však člověk na chvíli zamyslí, začne mu docházet, že monoidy najdou uplatnění na mnoha místech. A když už ne monoidy jako takové, tak aspoň jejich chudší příbuzní pologrupy (anglicky semigroup), které vypadají úplně stejně, jen jim chybí nulový prvek1.

Například operace sčítání s nulovým prvkem 0 tvoří monoid. To samé pro násobení a 1, spojování řetězců a prázdný string, spojování seznamů a prázdný seznam, nebo sjednocení množin a prázdnou množinu. Protože jde o monoidy, všechno se chová uniformě, regulárně a podle diktátu pravidel o asociativitě a identitě.

Několik jednoduchých monoidů:

typ plus zero
Int + 0
Int * 1
String + ""
List ++ List()
Set union Set()
Boolean and true
Boolean or false
Int max -Infinity
Int min +Infinity

To ale zdaleka není všechno. Monoidy je totiž možné kombinovat a komponovat a tvořit větší z několika menších. Jestliže nám například dva páry hodnot (ve Scale jde o Tuple2) a chci je sečíst, udělám to takhle:

(a₁, b₁) ⊕ (a₂, b₂) = (x, y)

Použiji příslušný monoid pro první komponenty a zvlášť monoid pro druhé komponenty.

x = a₁ ⊕ a₂
y = b₁ ⊕ b₂

Takže když páry mají typ (Int, String), budu sčítat první komponenty a spojovat druhé komponenty. O tom, jaký monoid se použije rozhoduje typ argumentů, které Algebird používá k výběru příslušné typeclass, která je předána jako implicitní argument.

val p1 = (1, "asd")
val p2 = (2, "xyz")

val (a, b) = Monoid.plus(p1, p2)

Kompilátor podlední řádek rekurzivně expanduje na něco ve stylu

Monoid.plus(p1, p2)(Monoid.monoid2(IntMonoid, StringMonoid))

Monoid existuje pro každý Tuple2, jehož obě komponenty mají monoidy. To samé pro Tuple3, Tuple4, atd. Stejně tak platí, že pokud mám jednoduché monoidy, můžu z nich konstruovat větší a složitější jako například Option, Either, Seq, Map nebo Function1.

Typeclass pro Option je definován takto:

class OptionMonoid[T](implicit semi: Semigroup[T]) extends Monoid[Option[T]] {
  def zero = None
  def plus(left: Option[T], right: Option[T]): Option[T] = {
    if (left.isEmpty) {
      right
    } else if (right.isEmpty) {
      left
    } else {
      Some(semi.plus(left.get, right.get))
    }
  }
}

Nulový prvek je None, plus sečte vnitřky Option pologrupou pro typ T, pokud jsou oba Some, jinak vrátí ten, který není None.

Podobně funguje Map: plus spojí dvě mapy a když obě obsahují určitý klíč, sečte jejich obsah příslušným monoidem.

Algebird se sekvencemi typu Seq zachází jako s tuple typy – nespojuje kolekce, ale sčítá odpovídající hodnoty. Sekvence typu List spojuje.

Jako příklad budu uvažovat, že mám kolekci čísel

val numbers: Seq[Int] = ???

Když ji chci sečíst, stačí zavolat:

val sum = Monoid.sum(numbers)

Výchozí monoid pro typ Int je sčítání + 0.

Pokud chci průměrnou hodnotu, musím obalit číslo typem AveragedValue. Algebird pro něj najde odpovídající typeclass, která implementuje plus jako průměrování Intů.

val average = Monoid.sum(
  numbers map { num => AveragedValue(num) }
)

Pokud chci zjistit více agregací, stačí mapovat výchozí data na n-tici požadovaných typů a Algebird udělá zbytek. Součet, minimum, maximum, průměr a množinu všech unikátních hodnot můžu v jedné iteraci zjistit například takto:

val (sum, min, max, avg, uniqueNumbers) =
  Monoid.sum(numbers map { x =>
    (x, Min(x), Max(x), AveragedValue(x), Set(x))
  })

Metoda sum(vs: TraversableOnce[T]): T definovaná na traitu a objektu Monoid sečte všechny prvky předané kolekce příslušným monoidem. V podstatě dělá jen values.foldLeft(zero)(plus).

Pochopitelně tohle není všechno. Monoidy se dají kombinovat mnohem vynalézavěji ve snaze spočítat komplikovanější agregace. Představte si, že mám access log webserveru a chci detekovat roboty. Začnu třeba tím, že spočítám kolik přístupů pochází z jednotlivých IP adres.

val records: Iterator[AccessLogRecord] = readAllLogRecords()

val map: Map[IpAddress, Int] =
  Monoid.sum(records map { r =>
    Map(r.clientIpAddress -> 1)
  })

Výsledkem je mapa, kde klíče jsou IP adresy a hodnotami je počet hitů.

Dobré. Ale co kdybych teď chtěl tyto počty pro každý měsíc zvlášť? Toho se dá docílit jednoduše:

val map: Map[Month, Map[IpAddress, Int]] =
  Monoid.sum(records map { r =>
    Map(yearAndMonth(r.date) -> Map(r.clientIpAddress -> 1))
  })

A co třeba statistiky o tom, kdy se daná IP adresa poprvé a naposledy objevila a kolik z ní bylo provedeno dotazů dohromady a v každou denní hodinu?

val res: Map[IpAddress, (Min[Date], Max[Date], Int, Map[Hour, Int])] =
  Monoid.sum(records map { r =>
    Map(r.clientIpAddress -> (Min(r.date), Max(r.date), 1, Map(hourOf(r.date) -> 1)))
  })

Monoidy neexistují jen pro jednoduché součty, minima, maxima a průměry, ale i pro komplikovanější objekty jako je například Bloom filter, HyperLogLog, Count-min sketch nebo top-k. Díky tomu můžu dělat komplikované agregace, které zahrnují pravděpodobnostní algoritmy, ale chovají se stále stejně jako jakýkoli jiný monoid a jsou stále komponovatelné.

Algebird nabízí tyto možnosti: DecayedValue, DecayedVector, Interval, Moments, BloomFilter, HLL, CMS, SpaceSaver, SketchMap, MinHash, HashingTrick, TopK.

Když budu chtít zjistit 100 dotazů, které vedly k největším odpovědím, můžu to udělat takhle:

// top-100 řazeno sestupně podle bytesSent
implicit val topk = new TopKMonoid(100)(Ordering.by[AccessLogRecord, Int](_.bytesSent).reverse)

val res: List[AccessLogRecord] =
  Monoid.sum(records map { r => topk.build(r) })

Nebo třeba odhadnout počet unikátních IP adres a unikátních uživatelů v každém měsíci pomocí HyperLogLogu, který potřebuje konstantní množství paměti nehledě na to, jaká je kardinalita množiny unikátních hodnot:

implicit val hll = new HyperLogLogMonoid(bits = 12)

val res: Map[Month, (Double, Double)] =
  Monoid.sum(records map { r =>
    Map(yearAndMonth(r.date) -> hll.toHLL(r.clientIpAddress), hll.toHLL(r.userId))
  }).mapValues { case (ips, users) => (ips.estimatedSize, users.estimatedSize) }

Jak je vidět, použití je vždycky stejné: nejdřív transformuji data do cílového tvaru a pak rozjedu sumaci a mechanismus monoidů podpořený implicitními parametry Scaly se postará o zbytek. Nezáleží přitom jak triviální nebo komplikovaná jsou data, když pro ně existuje monoid, dají se sečíst.


To, že dvojice nulového prvku a asociativní operace se vyskytuje často svědčí časté použití metod/funkcí jako je foldLeft:

xs.foldLeft(0)((a, b) => a + b)

V tomto případě nulový prvek a asociativní operaci monoidu. Jde tedy o:

xs.foldLeft(M.zero)(M.plus)

neboli:

Monoid.sum[M](xs)

Z rozhraní monoidu dále vyplývá, že je možné provádět výpočet dávkově, na (potenciálně nekonečném) proudu dat, paralelně a je možné dělat snapshoty. To všechno díky asociativitě binární operace.

Streamování je jasné – když dostanu nová data, přičtu je k současnému výsledku a všechno jede bez problémů. Musím jednom zaručit že nedošlo ke změně pořadí v jakém položky přicházejí, protože operace není nutně komutativní. Možnost snapshotů je pak duál této vlastnosti.

Paralelizace vyplývá z faktu, že můžu výraz libovolně uzávorkovat.

Například

a+b+c+d+e+f+g+h+i+j+k+l+m+n+o+p+q+r+s+t+u+v+w+x+y+z

můžu uzávorkovat takhle

(a+b+c+d+e+f)+(g+h+i+j+k+l+m)+(n+o+p+q+r+s+t)+(u+v+w+x+y+z)

každou závorku spočítat v jednom vlákně a pak sečíst mezivýsledky

AF + GM + NT + UZ

Opět je potřeba zaručit jenom, že nedojde ke změně pořadí položek a mezivýsledků (nebo zaručit komutativitu2) a paralelizace je zcela bezproblémová a koncepčně zadarmo.

Tady skončím, protože tohle stačí jako odpověď na tu nikdy nevyslovenou otázku z let dávno minulých. Je překvapivé, že tak málo může poskytnout tolik garancí a flexibility, a že tolik reálných a užitečných věcí, strukturou odpovídá něčemu, co na první pohled vypadá jako akademická podivnost popsaná několika řeckými písmeny a zcela nepřístupná reálnému světu pragmatiků, kteří se dřou pod praporem worse is better.

V podobném duchu nečekané abstrakce, která se objevuje na nečekaných místech se nese přednáška Regexes, Kleene algebras, and real ultimate power!. Člověk nemusí pochopit všechny detaily, aby cítil nepopiratelnost představované abstrakce.


Jako člověk, kterému se neštítí počítání taktů procesoru se musím zmínit ještě o jedné věci – výkonu. Algebird byl navržen s koncepční čistotou, která ne nutně musí vést k rychlému výslednému programu. Sumace pomocí monoidů uvažuje, že všechny objekty jsou neměnné a při každém součtu se vytvoří další, což vede k poměrně divokým alokacím. To je ještě horší problém, když sčítám mapy tuplů vnořené do jiných map tuplů. Při každé iteraci je velká část této struktury znovu vytvořena. Ne celá, protože persistentní datové struktury používají strukturální sdílení pod kapotou. Přesto jde o řádově víc alokací, než by dělal kód, který manipuluje měnitelnými strukturami. Výsledek může být až několikrát pomalejší. Ale na druhou stranu je sumace prostřednictvím monoidů nesrovnatelně jednodušší (hlavně, když se do hry dostanou hluboce vnořené struktury) a bezpečně paralelizovatelná. Situaci také pomáhá fakt, že ve scale jsou immutable mapy a množiny ručně specializovány pro 1 až 4 položky. Alokace malých map/množin je tedy relativně levná, protože není třeba vytvořit celé HAMT. Technicky by bylo možné provést stream fusion (například pomocí maker), které by odstranilo nutnost alokace dočasných a efemérních objektů.

Neříkám, že kvůli tomu jsou monoidy špatné nebo pomalé. Jde o jednu věc, o které je třeba se zmínit.


  1. Fakt, že pologrupy postrádají nulový prvek lehce komplikuje některé operace. Když chci sečíst libovolnou sekvenci dat, pologrupa si s tím neumí poradit, protože libovolná sekvence může být prázdn. V takovém případě by se vyplatilo vrátit nulový prvek. V Algebirdu je proto na třídě Semigroup deklarována metoda sumOption(xs: TraversableOnce[T]): Option[T] a sum zcela chybí.
  2. Mnoho operací je komutativních, protože jsou zároveň idempotentní: min, max, sjednocení množin, logické operace, konstrukce Bloom filteru a HyperLogLogu. Jiné nejsou idempotentní, ale přesto komutují: sčítání, konstrukce Count-min sketch.

Dále k tématu:

Flattr this!

Posted in Funkcionální programování, Scala | 1 Comment

Von Neumannovy lži

Asi takhle: Existují dvě široké rodiny procesorů. Na jedné straně Von Neumannovy stroje a na druhé data-flow mašiny. První jmenovaná zpracovává jednu instrukci za druhou, druhá k programu přistupuje jako ke grafu závislostí, kdy každá operace je závislá na některých předcházejících a jakmile jsou všechny závislosti splněny, může být operace provedena.

Von Neumannova architektura je na papíře horší a mnohem méně ambiciózní než data-flow, protože je z povahy sekvenční, kdežto data-flow architektury jsou přirozeně paralelní. Jak dokládá historie, nakonec jsme skončili ve spárech Von Neumanna i přesto, že data-flow nebyl jen divoký sen poblázněných akademiků a tyto stroje se skutečně vyráběly. Na první pohled to může vypadat jako další triumf worse is better a vítězství průměrnosti. Pravda je však komplikovanější a má v sobě něco málo dramatické ironie.

To by jako úvod stačilo, teď je čas se po hlavě vrhnout do zákopů.

Budu uvažovat, že jsem napsal následující funkci (jde o kus kódu z mých experimentů s řazením):

def packedTerminatedByteSlice(arr: Array[Byte], from: Int, to: Int): Long = {
  var i = from
  var j = 0
  val end = math.min(arr.length, to)
  var res = 0L

  while (i < end) {
    res |= ((1 << 8) | arr(i)) << j
    i += 1
    j += 9
  }

  res
}

Jak by tato funkce byla vykonána na hypotetické data-flow mašině? Zajímá mě jen smyčka na konci funkce.

Je vidět, že v těle smyčky existují tři nezávislé cesty, které může data-flow stroj provádět najednou. Zatímco dělá logické ORy a SHIFTy, aby aktualizoval proměnnou res, může paralelně inkrementovat i a j. Napříč iteracemi se otevírají další možnosti: Podmínka smyčky i < end závisí jen na proměnné i, je možné ji tedy testovat, zatímco stroj provádí operace minulé iterace. Ze schématu je vidět, že jeden běh smyčky zabere tolik času, kolik potřebuje nejdelší cesta v grafu (modulo implementační detaily).

Von Neumannova architektura je proti tomu mnohem jednodušší. Program je přeložen na lineární sekvenci instrukcí, které jsou prováděny v řadě jedna za druhou.

Moje funkce by se mohla přeložit například na následující pseudo-assembly:

cond = cmp i end
jump_if_less_then cond
tmp1 = mov arr(i)
tmp2 = or 256 tmp1
tmp3 = shift tmp2 j
res  = or res tmp3
i    = add i 1
j    = add j 9
jump back

Kanonický Von Neumann načte jednu instrukci, provede ji a přesune se na další. Tělo smyčky obsahuje 9 instrukcí a to by mělo určit dobu jejího běhu1.

Návrháři architektur jsou však vynalézaví lidé, věčně nespokojení s tím, co současný hardware dokáže, a neustále přemýšlejí, jak starého Von Neumanna vylepšit. Došlo ke zvyšování taktů, což vedlo nejdříve k několika-taktovým instrukcím a pak k pipelinování. Ale stále šlo o stejnou architekturu – jedna instrukce v jeden čas, i když se různé fáze různých instrukcí překrývaly. Změnila se jediná věc: Nejednou bylo třeba sledovat závislosti pipelinovaných instrukcí. Pokud operace hned požadovala výsledek té předchozí, musela čekat, než ta předchozí dokončí svojí práci a zapíše výsledek do registru2. Kdyby nečekala, přečetla by nevalidní data. Hardware musel na pár taktů zastavit práci a v pipeline se vyskytla bublina (pipeline bubble).

a = x + 1
b = a + 1 // závisí na a

a: [ fetch | decode | read operands | execute | write result ]
b:         [ fetch  | decode        | ------- | ------------ | read operands | execute | write result ]

Pozorný čtenář si jistě všiml, že poslední dva součty ve výpisu instrukcí ukázkové funkce jsou na sobě zcela nezávislé. Je tedy možné je odpálit najednou bez komplikované data-flow mašinérie, která sleduje závislosti celého programu. Stačí jen trochu křemíku, který analyzuje vztahy sousedících operací. Tohle zjištění vedlo ke vzestupu superskalárních/wi­de-issue strojů – stále Von Neumannovských architektur, jen rozšířených o schopnost najednou odbavit nezávislé sousedící instrukce, pokud hardware má dostatek prostředků (v mém případě dvě ALU). Může jít o dynamickou vlastnost, kdy třeba X86 procesor analyzuje proud instrukcí za běhu a nebo statickou jako v případě VLIW nebo Itania.

Wide-issue jak jsem ho popsal využívá paralelismu jen v omezených případech. Nemůže paralelně provádět nezávislé instrukce mimo pořadí programu. Například tento kousek kódu

tmp1 = add input 1
x    = add tmp1 1
y    = add input 2

má vzájemně nezávislé operace na prvním a třetím řádku, ale druhý řádek závisí na výsledku toho prvního. Aby hardware mohl paralelně provést nezávislé instrukce, musí dělat věci mimo pořadí programu, tedy out-of-order (OOO). Tímto směrem směřovala evoluce ctihodného Von Neumannova stroje. OOO hardware sleduje vzájemné závislosti poměrně velkého okna instrukcí3 a snaží se naplno vytížit všechny dostupné funkční jednotky pomocí spekulativní exekuce. Udržuje si při tom dostatečné množství stavových informací, aby se dokázal vrátit zpátky do konzistentního stavu, kdyby se spekulace ukázaly jako příliš optimistické (např. kvůli špatné branch prediction).4

Všimli jste si toho taky? Out-of-order procesor vypadá skoro jako data-flow stroj.

Ironií osudu je, že postupnou evolucí se Von Neumann postupně přibližoval data-flow architekturám a došlo to tak daleko, že interní chování OOO stroje je téměř k nerozeznání od jeho data-flow bratrance. Rozdíl je jen v tom, že OOO stroj vytváří graf závislostí za běhu na poměrně omezeném okně instrukcí, zatímco data-flow to dělá na celém programu v době kompilace.

OOO stroje do segmentu osobních počítačů vpadly spolu s Pentiem Pro v roce 1995. Von Neumanovu architekturu tedy téměř nikdo z nás nepoužívá už dvacet let.


Dále k tématu:

  1. Pokud tedy budu uvažovat jednoduchý model, kdy každá instrukce trvá stejně dlouho a odmyslím si všechny nežádoucí efekty jako jsou rozdíly mezi rychlostí cache a RAM.
  2. Tyto případy částečně řeší tzv. Operand forwarding, který kromě zápisu do příslušného registru, data přesměruje přímo do následující instrukce.
  3. Například na Haswelech je to 192. Skylake generace Intelích CPU to rozšířila na 224. Nejde o x86 instrukce, ale o dekódované mikro-operace.
  4. Způsoby zotavení ze špatných spekulací jsou popsány například ve Fast Branch Misprediction Recovery in Out-of-order Superscalar Processors

Pozn:

Flattr this!

Posted in CPU | 5 Comments