Čím více se věci mění, tím více zůstávají stejné

Evgeny Morozov v knize To Save Everything Click Here mluvil o epochalismu a ahistorickém myšlení. Epochalismus říká, že žijeme ve výjimečných časech a současná doba je diametrálně odlišná od všech minulých a nic staré už neplatí. Ahistorické myšlení s epochalismem velice úzce souvisí – v podstatě jde o systematické ignorování historie a všeho, co jsme se z ní mohli naučit.

Přestože Morozov psal o „Internetu“ a digitální debatě obecně, stejný trend lze vycítit v IT komunitě. Zdá se, že dominantní myšlenkový proud tvrdí, že se všechno neustále mění a vyvíjí, že co jsme se naučili včera bude zítra zastaralé, že co jsme věděli přestane platit a všechno se musí co nevidět od základů změnit. Cynik by ukázal na seznam Javascriptovch frameworků a namítl, že nic jako pokrok neexistuje a že jen opakujeme stejné chyby a občasný skok vpřed nastane, když znovuobjevíme zapadlou akademickou práci ze sedmdesátých let.

Já tak daleko nechci zajít, protože popírat samotnou myšlenku pokroku je šílené. Zdá se však, že iluze totální disrupce je v oboru, o kterém by cynik prohlásil, že je poháněn vpřed výhradně nesmysly, marketingem a popíráním minulosti, nevyhnutelná.

S tímhle na mysli jsem se začal zajímat, kdy a kde byly vynalezeny všechny techniky organizace a architektury používané v současných procesorech. Kdyby ahistorická hypotéza platila, všechno na čem jsou soudobá CPU postavená, se diametrálně liší od všeho minulého a jde jen o dočasnou bublinu, která co nevidět splaskne. Pokud je to pravda, nemá valný smysl současný hardware sledovat.


Všechno začalo dvěma systémy: ILLIAC II (1962) a IBM Stretch (1956–1961), jejichž tvůrci přišli s naprostou většinu architekturálních technik a způsobů organizace procesorů.

Tyto dva stroje uvedly cache, instrukční pipeline, multitasking, virtuální paměť, ochranu paměti, osmibitový bajt, prefetch a prokládání paměti (který je bohatě používán např. v GPU).

Po těchto naprosto zásadních strojích následoval systém CDC 6600 (1965) a výzkumný projekt IBM ACS-1 (později přejmenovaný na ACS-360, 1961–1969). Ty uvedly všechno ostatní. Jmenovitě: dekódování a vykonání několika instrukcí v jednom taktu (superskalární CPU), branch target buffer, multithreading implementovaný v hardware, out-of-order execution, přejmenovávání registrů, predication.

Současné procesory obsahují tohle všechno a nic víc. Poslední generace hardware jsou jen sny designérů šedesátých let dovedené k dokonalosti.

Je však třeba dodat, že i když byla nějaká technika přivedena k životu před padesáti lety, neznamená to, že byla celou dobu aktivně využívána. Velká část postupů byla na dlouhá léta opuštěna, protože se dařilo výkon dostatečně zvyšovat prostým nárůstem frekvence. Nebylo třeba dělat nic chytrého, stačilo se nechat unášet vlnami Moorova zákonu a stále pokročilejší litografií. Až s postupem času (ale stále dlouho před zastavením růstu frekvence1) začali návrháři čipů využívat rostoucí množství tranzistorů k agresivním implementacím chytrých a pokročilých technik.

Ve zbývajících odstavcích projdu jednotlivé architektonické vzory a jejich stručnou historii.

Instrukční pipeline

Instrukční pipeline byla použita už ve strojích ILLIAC II a IBM Stretch. První nasazení je však o mnoho let starší a sahá až k mechanickým a elektromechanickým počítačům Z1 (1939) a Z3 (1941). Pipeline poté začala být používána v superpočítačích sedmdesátých let a je spjata s přechodem od jedntaktových instrukcí k mnohataktovým instrukcím.

Cache

Už ILLIAC II měl rychlou a pomalou paměť. Pomalá paměť měla přístupovou dobu 2 µs, rychlá, které by se dnes říkalo cache, 0.25 µs.

Rozdíl mezi rychlostí procesorů a pamětí začal narůstat v osmdesátých letech. Paměť, která byla předtím dostupná v jednom taktu, přestávala stíhat CPU a čtení a zapis začal trvat víc jak jeden takt.

První cache byly TLB, které se objevily krátce po příchodu virtuální paměti ve strojích GE 645 (1968) a IBM 360/67 (1967).

Datová cache jako taková se poprvé objevila ve stroji IBM S/360 M85 v roce 1968.

Dedikovaná instrukční cache měla premiéru v Motorole 68010 (1982), do které byl přidán „loop mode“ pro 2 instrukce těsné smyčky. Pozdější Motorola 68020 (1984) přidala 256 bajtovou I-cache, Motorola 68030 (1987) přidala stejně velkou D-cache. První x86 procesor s instrukční cache bylo dvacetimegahertzové 386 (1987).

Superskalární CPU

Už CDC 6600 z roku 1965 byl superskalární stroj. Prvními komerčními jednočipovými superskalárními procesory byly až Motorola MC88100 (1988), Intel i960 (1989) a AMD 29000 (1990, viz). Pentium se v roce 1993 stalo prvním superskalárním x86 procesorem. Více méně všechny procesory vyrobené přibližně po roce 1998 jsou superskalární.

Nx586, Pentium Pro a AMD K5 byly první x86 superskalární procesory, které překládaly složité CISC instrukce na jednodušší RISC operace, kterých pak mohly vykonat víc najednou.

Out-of-order

Prvním OOO strojem byl CDC 6600 z roku 1965. Avšak podle dnešních standardů nešlo o plnohodnotné OOO, protože používal scoreboarding, který dovoloval pouze dokončení (retire) instrukcí out-of-order, začínat musely stále v pořadí programu. To umožňovalo začít rychlé a pomalé operace najednou s tím, že ty pomalé nevytvořily bublinu v pipeline, ale hardware je prováděl souběžně a potom přerovnal jejich výsledky tak, aby odpovídaly sekvenčnímu programu.

Skutečný OOO stroj uvedl až big blue a byl jím IBM 390/91 z roku 1966 implementující Tomasulo algoritmus, který sleduje závislosti mezi instrukcemi a dynamicky je plánuje mimo pořadí. Ale i tady to mělo jeden háček: OOO logika byla použita jen pro floating point část procesoru. To bylo nejspíš proto, že FP operace byly drahé a FP logika zabírala velkou část čipu a bylo třeba zařídit, aby byla za každou cenu efektivně využita.

Pak se třicet let nic nedělo a prvním moderním OOO mikroprocesorem se stal až POWER1 (1990), který také OOO aplikoval jen na floating point operace. O tři roky později IBM uvedl PowerPC (1993), který už uměl plné OOO. To rozpoutalo bouři a mnoho architektur začalo přecházet: Sparc (1995), Pentium Pro (1995), R10000 od MIPS/SGI (1996) PA-RISCový PA-8000 od HP (1996), AMD K5 (1996) a nakonec Alpha (1998).

Branch prediction

Už IBM Stretch měl jednoduchý statický branch predictor, který vždycky předpokládal, že podmíněný skok nebude použit a kód propadne (predict untaken)2. Po tomto experimentu IBM nepoužívala branch predictor ve velkých systémech až do roku 1985. Pravděpodobně to nebylo třeba, protože cena špatného odhadu byla s krátkou pipeline malá nebo byly k dispozici predikované instrukce, případě delay sloty

Branch prediction se objevuje v systémech jiných společností: Burroughs B4900 (1982) nebo VAX 9000 (1989). První komerční RISCy (1986) po vzoru Stretche a měly statický prediktor predict untaken. Branch prediktor začal být velice důležitý po roce 1993 pro superskalární procesory s hlubokou pipeline, kde cena špatného rozhodnutí začínala narůstat. Pentium, Alpha 21064, R8000 a POWER také začaly nasazovat branch predictior.

Multicore

Překvapivé je, že i několikajádrové procesory nejsou tak nové, jak by se zdálo. První kousek vyrobil v polovině osmdesátých let Rockwell, když zkombinoval dva 6502 čipy do jednoho pouzdra. Po roce 2000 se staly mnohojádrové procesory běžnými artikly.

SMT/HyperThreading

Hardwarový multi-threading (SMT) byl poprvé implementován v prototypu ACS-360 (1968), který nebyl nikdy dokončen.

Procesor Alpha 21464 (ohlášený 1999, plánovaný na 2003, nakonec zrušen), měl být prvním SMT mikroprocesorem, měl zvládat 8-wide issue (8 instrukcí v jednom taktu) a 4 souběžná vlákna v každém jádře. Pentium 4 (2002) se stalo prvním moderním desktopovým CPU s podporou SMT. SMT v podání Intelu (překřtěné na HyperThreading) zvládá jen dvě vlákna na jednom jádře. Další na řadě byl POWER5 (2004) od IBM.

SIMD

Vývoj vektorových procesorů, které spadají do kategorie SIMD, začal v šedesátých letech (např. ILLIAC IV). Konkrétní stroje se začaly objevovat v letech sedmdesátých a šlo především o různé architektury Cray (viz).


Jak je vidět, všechny techniky organizace a architektury procesorů jsou velice staré. Naprostá většina zažila debut na přelomu padesátých a šedesátých let. Cynik by mohl namítat, že se v tomto oboru od osmdesátých let nestalo nic skutečně nového a nebyl by daleko od pravdy. Dalo by se říct, že v procesorech je za posledních padesát let to samé, jen je tam toho víc.6 Protože jsou tyto principy s námi od samého počátku, dá se předpokládat, že na nich něco bude a budou mít určitou univerzální platnost. Jinými slovy: když budeme optimalizovat pro současný hardware, neznamená to nutně, že podléháme diktátu podivností poslední várky křemíku z továren Intelu, ale široké množině procesorů, které mají velkou cache, hlubokou pipeline, chytrý branch predictor, agresivní OOO engine a neuvěřitelně pomalou DRAM. Do této kategorie spadá skoro všechno od největších serverových bestií až po mobilní procesory. Ona poslední generace mobilních procesorů se skoro vůbec neliší od těch desktopových a třeba papírová specifikace Apple Cyclone se až nepříjemně podobá číslům Haswellu.3

Některé z těchto principů nejsou jen výsledky snah hardwarových designérů napravit hříchy našich otců4, ale důsledky fyzikálních principů. Cache je například důsledkem toho, že v prostoru může být jen omezené množství věcí blízko procesoru (rychlá cache) a obrovské množství daleko (pomalá DRAM). Datová závislost je potom manifestací kauzality, kdy pointer neukazuje jen na místo v paměti, ale zároveň udává pořadí operací.

What's the score here? What's next?

Co by se muselo změnit, aby staré pořádky přestaly platit?

Na paměti s malou latencí nejspíš můžeme zapomenout. Tahle bitva je už dávno prohraná. Místo toho se můžeme dočkat strojů s velikou paměťovou propustností. Každý jeden přenos zabere dlouhou dobu, ale hardware jich může provádět hodně najednou, podobně jako to dělají grafické karty nebo ThreadStorm architektura. Tenhle datově paralelní model funguje nejlépe, když jsou operace relativně přímočaré a vzájemně nezávislé. To ale není pravda o obecném kódu, který je plný podmíněných skoků, polymorfismu a často vyžaduje jemnou synchronizaci a koordinaci. GPU mají problém s podmíněnými skoky, alokacemi a rekurzí, tedy s tím, čeho je general purpose kód plný. Realita bude tedy nejspíš taková, že obecné procesory budou mít u sebe paměť s velkou propustností, ale jednotlivá jádra budou celkem obyčejná jádra podobná těm, jaké máme v současnosti, jen možná o něco jednodušší, možná s širším SMT nebo hardwarovým přepínáním vláken.

Pokrok můžou přinést levné nevolatilní paměti (NVRAM), které fungují jako RAM, ale data přežijí výpadek napájení a restart systému (jako 3D XPoint od Micronu/Intelu nebo memristorová paměť od HP5). Mnoho lidí je z nich nadšených, protože NVRAM můžou změnit úplně všechno (co se týká RAM). Na druhou stranu také nemusí změnit vůbec nic. To, že data přežijí pád sytému, znamená také, že porušená data přežijí restart a chyby můžou přetrvat. Někdo může jásat nad tím, že nebude třeba dělat nic speciálního pro uložení souboru, protože data jsou v paměti a paměť je persistentní. Pokud tohle přijmu, musím uvažovat, co se stane, když systém spadne v průběhu jedné operace. Najednou můžou být data zapsána částečně a kdo ví, co to bude dělat, až systém znovu naběhne. V tomto případě je třeba mít transakční protokol úplně všude a ne jen v místech, které se starají o ukládání dat. Z tohoto důvodu možná nějaká paměť nebude tolik persistentní jako jiná. Když přihlédnu k prezentovaným číslům, která naznačují, že NVRAM může být rychlá asi jako DRAM, ale přesto o něco pomalejší, některé systémy budou mít na palubě jak NVRAM tak i DRAM. Nakonec to nemusí být technologie, která nás všechny spasí, ale jen jeden druh z celého spektra pamětí: páska, HDD, SDD, NVRAM, DRAM, SRAM.

Když zajdu do jiného oboru a podívám se na řadící algoritmy je situace až nepříjemně podobná: Všichni stále používáme quicksort, na který v roce 1959 přišel Tony Hoare (publikováno v roce 1961 jako Algorithm 64: Quicksort), merge sort, jehož autorem byl sám John von Neumann (1945), nebo Radix sort, jehož kořeny sahají až do roku 1897. Kromě nich nepoužíváme skoro nic jiného. To ale neznamená, že se nic neděje. Všechny zmíněná algoritmy jsou postupně vylepšovány, vylaďovány a přizpůsobovány hardwarové realitě doby, jako využití SIMD bitonických sítí pro merge sort, skosený quicksort nebo různá schémata paralelizace všech sortů. Ale pořád platí, že znalost základních verzí těchto algoritmů nikdy nezastarala.

Naproti tomu svět databází se otřásá v základech. Tedy to tak aspoň vypadá na první pohled. Současné RDMS v některých případech přestávají stačit požadavkům na ně kladeným a proto se začaly ozývat hlasy k zapovězení relačních databází a přechodu na něco drasticky jiného a lepšího. Dokonce i Stonebraker tvrdí, že doba jedné databáze pro všechno je u konce a je efektivnější přejít na in-memory databázi na OLTP, sloupcovou na OLTP a něco jako Hadoop/Spark na ostatní věci. V tomto kontextu není těžké pochopit hlasy, které tvrdí, že relační databáze jsou zastaralé a morálně špatné. To ale nic nemění na tom, že byly více než dostačující posledních 40 let, od doby, kdy Codd publikoval seminální paper o relační algebře. Navíc techniky použité pod kapotou se příliš neliší mezi (některými) relačními dinosaury a (některými) moderními key-value NoSQL datastory. Konec konců chci někam uložit hromadu bitů a chci ji taky dostat rychle zpátky.

V dlouhodobém měřítku všechno jednou zastará a jediné, co můžeme chtít, je pár let jízdy plnou rychlostí, než to strhneme do škarpy a zmizíme v explozi slávy. Proto ale nemá smysl tvrdit, že dočasná znalost nemá žádnou hodnotu. Jestliže to, co se dneska naučíme o hardware, bude platit následující dekádu, bude to víc než bychom si mohli přát. Onen pomyslný cynik by jistě ukázal na hořící hromadu Javascriptových frameworků a dodal by, že je to víc než bychom si zasloužili.


Dále k tématu:

Pozn:

  1. viz známý článek The Free Lunch Is Over.
  2. Z toho se dá usoudit, že předpokládal, že smyčky budou začínat skokem ven a na konci bude nepodmíněný skok na začátek smyčky. Kdyby předpokládal opak (predict taken), znamenalo by to, že smyčce bude končit podmíněným skokem na začátek.
  3. O tom také svědčí fakt, že mnoho firem plánuje uvést vlastní ARM serverové čipy. Jsou to třeba AMD, Qualcomm, Tilera, EZchip, Cavium, Aplied Micro nebo Broadcom.
  4. viz bizarní kódování x86, které vede k masivním CISC → RICS dekodérům, L0 cache, LSB a dalším snahám jak rozlámat instruction stream a nakrmit hladové funkční jednotky.
  5. Pro úplnost se sluší dodat, že NVRAM je dostupná už teď, jen není vůbec levná (MRAM, FRAM). Příslib NVRAM je v tom, že bude jak nevolatilní, tak i velká.
  6. Podle přednášky The Chip Design Game at the End of Moore's Law se mezi roky 1980 a 2010 zvýšil takt procesorů 3500× a změny v architektuře a mikro-architektuře mohly přidat další padesátinásobné zvýšení rychlosti. Z toho plyne, že se skutečně něco děje a všechny ty staré principy se skutečně vylepšují a vylaďují k dokonalosti.

Flattr this!

Posted in CPU, DB | Leave a comment

How to sort out your life in O(n) time

slides

Flattr this!

Posted in Uncategorized | Leave a comment

Jak řadit v lineárním čase, křísit mrtvé a dosáhnout osvícení

Když se řekne řazení, většina lidí, kteří napsali víc než jednu smyčku a dva ify se zeptá: „A je to nutné?“ Řazení nemá pověst nejrychlejší operace pod sluncem a je preferováno, když se tomu dá vyhnout např. hashováním. Někdy je to ale nutné a v tom případ je třeba zatnout zuby a připravit se na nějaký ten rozpálený křemík.

Jedinci, kterým na zdi visí portrét Donalda Knutha si okamžitě vybaví, že řazení má v nejlepším případě časovou složitost O(n log(n)) a vzpomenou si na všechny ty quicksorty, merge sorty, heap sorty a různé monstrózní hybridy.

Pravda je ale o něco komplikovanější: Řadící algoritmy založené na porovnání potřebují provést řádově n * log(n) porovnání. Existují však i jiné způsoby řazení, které jsou založené na frekvenci, rozložení a histogramech, a ty nepotřebují porovnávat nic s ničím. Stačí jim lineární čas, a co víc: Jsou zatraceně rychlé.

Jedním z nejjednodušších je counting sort, který funguje tak, že vstupní pole jednou projdu, zaznamenávám kolikrát jsem na který element narazil a na konci položky vyplivnu ve požadovaném pořadí.

Může to vypadat například takhle (psáno v jazyce fast-and-loose Scala):

def countingSort(arrayToBeSorted: Array[Int]): Array[Int] = {
  val buckets = new ArrayBuffer[Int]()

  for (el <- arrayToBeSorted) {
    buckets(el) += 1 // increment counter
  }

  val result = new ArrayBuffer[Int]()

  for {
    i <- 0 until buckets.size // element
    _ <- 0 until buckets(i)   // count
  } result += i // append element

  result.toArray
}

Bystří si jistě už všimli, že slibovaná lineární složitost má pár nedostatků. Jak je vidět, je třeba mít po ruce pole, které je velké jako největší element ve vstupní kolekci. Co se stane, když bych se pokusil seřadit čísla 1 a 1844674407370­9551616? Musel bych přikoupit hodně paměti.

Counting sort má lineární časovou složitost, ale potřebuje prostor přímo úměrný rozpětí možných hodnot. To je ideální, když řadím například pole bajtů, ale v ostatních případech to má k ideálu velice daleko.

Řešením je radix sort. Ten funguje velice podobně, ale k řazení nepoužívá celou hodnotu najednou (tj. všechny dostupné bity), ale vždy jen její část. Tak například, když bych chtěl seřadit pole čtyřbajtových intů, radix sort může v první iteraci data rozdělit podle nejdůležitějšího bajtu. Vznikne tak 256 skupin (partition), kdy každá z nich obsahuje jen elementy shodující se v nejvyšším bajtu. V další iteraci provedu stejný krok pro každou skupinu zvlášť s druhým nejdůležitějším bajtem. Vznikne tak 256×256 podskupin, z nichž každá sdílí horních 16 bitů. Když udělám 4 rekurzivní iterace a spojím vzniklé podoblasti, dostanu seřazenou sekvenci. Stačí mi k tomu jen 4 iterace vstupním polem a pokud se nepletu, 4 lze považovat za konstantu.

Časová složitost radix sortu je O(kn), kde n je délka vstupu a k je velikost jedné položky. Když řadím pole primitivních typů nebo malých structů, k je malé a všechno běží velice rychle. Když se však rozhodnu radixovat dlouhé stringy nebo velké a komplikované (nedej bože nekonečné) objekty, k začne dominovat a situace je o poznaní méně růžová. Zvlášť s přihlédnutím k faktu, že log2(232), což je maximální rozsah 32 bitového intu, je jen 32, k tedy nemá moc místa pro manévrování. Stačí když radixovaný objekt trochu nakyne a najednou n log(n) představuje správnou cestu vpřed.


Majíce tohle všechno na paměti jsem napsal vlastní implementaci radix sortu ve Scale, abych na vlastní oči viděl jak rychlé je rychlé řazení. Můj kód je založený na moudrosti tohoto článku a jde o LSD radix sort. Neřadí tedy od nejvýznamnějšího bajtu, ale od toho nejméně významného. To se může zdát neintuitivní, ale funguje to, protože radix sort je stabilní algoritmus – více detailů ve zmiňovaném článku. Pro čtyřbajtové inty potřebuje provést 5 iterací nad daty – v jedné spočítá histogram frekvencí jednotlivých bajtů a ve čtyřech zbývajících přesouvá položky podle histogramů.

To však není jediný způsob. Existuje spousta dalších odrůd radix sortu: je možné řadit od nejdůležitějších bitů, od těch nejméně důležitých, může být in-place, může potřebovat extra O(n) místa, může být hybridní (použije se jiný algoritmus, když je jasné, že se nevyplatí radixovat), může být paralelní, může brát v potaz cache, může používat software managed buffers, atd. atd.

Výsledek? Jak řekl klasik: ho-ly-shit. Přesvědčte se sami:

Graf ukazuje rychlost řazení polí obsahujících náhodné inty mým radix sortem a funkcí java.util.Arra­ys.sort, což je hybrid mezi dual-pivot quicksortem, mergesortem a insertion sortem. Jak je vidět, radix sort běží v lineárním čase a navíc je velice rychlý. A to jde o neaktuální graf, který měří rychlost verze bez několika optimalizací3, a navíc do průměru započítává i počáteční pomalé běhy, kdy JVM nebyla zahřátá. Malé zpomalení s rostoucí velikostí dat není způsobeno asymptotami algoritmu, ale hardwarovými efekty. Všechny operace jsou rychlé, čtení je vždycky sekvenční, zápis je skoro sekvenční do 256 míst v paměti, které se pohodlně vejdou do cache5, neobsahuje žádné nepředvídatelné skoky, které by vedly k brach-miss – situace je zcela odlišná od quicksortu, který neobsahuje nic jiného než nepředvídatelné skoky.1

1GB dat (250M zcela náhodných intů) je na mém stroji možné seřadit za 4 vteřiny, a 1GB floatů za 2.9 vteřiny2. To odpovídá 14.9 ns, resp 10.8 ns na každý element pole. Když se celé pole vejde do cache, řazení je o něco rychlejší. Každému intu i floatu stačí 9.7 ns, což odpovídá rychlosti řazení přes 100M/s v jednom vlákně. Radix sort navíc umí poznat, když není třeba řadit podle určitého bajtu a přeskočí celý jeden řadící krok, což vede ke značnému zrychlení. Například když řadím pole obsahující jen čísla mezi 0 a 64k, dvě iterace z pěti jsou přeskočeny a je třeba jen 8 ns místo 14.9 ns na element.

Řazení floatů radix sortem je poměrně delikátní záležitost, která se spoléhá na to, že čísla mají bitovou reprezentaci podle specifikace IEEE754. Nejjednodušší možností je data seřadit jako integery a pak opravit pořadí záporných hodnot. O něco rafinovanější je floaty na vstupu transformovat, tak aby by je bylo možné řadit jako 32-bitová čísla bez znaménka a na konci je změnit zpátky na float. Nejlepší řešení je však tomu předejít a upravit offsety záporných radixů.


Ok, takže když je radix sort tak rychlý, k čemu se dá použít? K řazení, pochopitelně. Ale najde i uplatnění tam, kde se řazení zdá, jako jít s kanónem na vrabce.

Jedno takové místo je seskupování.

Když bych chtěl seskupit moře hodnot podle jejich klíče z páru klíč-hodnota (kde klíč i hodnota jsou čtyřbajtové inty zabalené do osmibajtového longu), mohl bych to s do očí bijící neefektivitou napsat v jednom řádku Scaly.

keyvals.groupBy(getKey).map { case (k, vs) => vs.map(getKey) }

To je pěkné, ale nešlo by to rychleji, nejlépe s minimem alokací? Šlo. S pomocí mapy, která obsahuje měnitelné buffery.

val buffers = mutable.Map[Int, mutable.ArrayBuilder.ofInt]()

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  buffers.getOrElseUpdate(k, new mutable.ArrayBuilder.ofInt) += v
}

buffers.toSeq.sortBy(_._1).values

Když zanedbám režii mapy (která může být zanedbatelná, ale také nemusí), za pozornost stojí poslední řádek: sortBy. Najednou jsem tam, kde jsem začal. Když je velikost každé skupiny malá, je jejich počet úměrný velikosti vstupních dat a konečné řazení bude velice drané.

Když vím, že klíče náleží do určitého rozpětí, můžou použít pole bufferů.

val buffers = new Array[mutable.ArrayBuilder.ofInt](length)

for (kv <- keyvals) {
  val k = getKey(kv)
  val v = getVal(kv)
  if (buffers(k) == null) {
    buffers(k) = new mutable.ArrayBuilder.ofInt
  }

  buffers(k) += v
}

buffers.filter(_ != null).map(_.result)

Jak je vidět, na konci není třeba nic řadit, protože výsledek je přirozeně seřazen podle klíče ve stylu counting sortu. Pokud mě trápí opakované alokace a kopírování dat v rostoucích bufferech a můžu si dovolit dva průchody daty, je možné v jenom průchodu spočítat velikosti jednotlivých skupin, alokovat pro ně pole velká tak jak je potřeba, a v druhém průchodu je naplnit. To může být v některých případech rychlejší a v jiných, kdy nechci přijít s OutOfMemoryError po funuse, nezbytné.

Tohle je lepší, ale nešlo by to… však víte… rychleji?

Tady do hry vstupuje radix sort. Co kdybych pole seřadil podle prvních 4 bajtů (tj. té části, která obsahuje klíč) a pak výsledek v jedné iteraci rozřezal na kousíčky. Nějak takhle:

val scratch = new Array[Long](keyvals.length)
val (sorted, _) = RadixSort.sort(keyvals, scratch, 0, keyvals.length, 4, 8, false)

var i = 0
(0 until numberOfKeys iterator) map { key =>
  var start = i

  // find end of the current group
  while (i < keyvals.length && getKey(sorted(i)) == key) { i += 1 }
  val end = i

  // empty group
  if (start == end) {
    null

  // copy values  into a separate array
  } else {
    val res = new Array[Int](end - start)
    var j = 0
    while (start < end) {
      res(j) = getVal(sorted(start))
      start += 1
      j += 1
    }

    (key, res)
  }
} filter (_ != null )

Na rozdíl od předchozích případů je kód o něco složitější, na druhou stranu nejde o pseudokód, ale skutečně fungující fragment. A co víc: Je skutečně rychlejší než všechny předchozí ukázky.

Proč?

Jedno vysvětlení je nasnadě: Pokud je skupin hodně a data jsou více méně náhodná, v předchozím případě každý nový klíč může vést ke několika drahým cache-miss. Jeden načte více méně náhodný pointer na buffer a druhý načte objekt bufferu a třetí konec pole v bufferu na které se bude vkládat. A už jsem zmiňoval, že jít do DRAM je drahé? Myslím, že ano. Radix sort tyhle problémy nemá, párkrát projede data, seřadí je a výsledek se nakonec ještě jednou proletí a rozseká na kousíčky. A všechno udělá velice rychle.


Dále k tématu:


  1. Špatně odhadnuté skoky představují tak velký problém, že někdy je quicksort rychlejší, pokud je záměrně zvolen špatný pivot. Když pivot data rozdělí fifty-fifty, je skok nepředvídatelný. Když je ale rozdělí třeba 25–75, procesor bude odhadovat, že skok půjde doprava a většinou se trefí. viz: An Experimental Study of Sorting and Branch Prediction. V tomto kontextu za pozornost také stojí: Branch Prediction and the Performance of Interpreters – Don't Trust Folklore a Branchless code sequences.
  2. Předpokládám, že v případě floatů JVM JIT zkompiluje kód s využitím SIMD instrukcí a proto je výsledek rychlejší. Nebo také může jít o korelace mezi bitovými vzory v reprezentaci floatů. Podrobnosti jsem nezkoumal.
  3. Jedna z optimalizací, která může za velice specifických situací vést k dvojnásobnému zrychlení je zpětná iterace jedné smyčky. V běžných situacích, kdy jsou vstupní data větší než cache přidá k dobru ~10%. Když je vstupní pole velké, nemá žádný vliv. Další optimalizaci, kterou jsem zkoušel, byly software managed buffers zmiňované zde a zde. Ty fungují tak, že data nezapisuji přímo do paměti, ale do malého bufferu, který se vejde do cache a do paměti ho zkopíruji jen, když se naplní. Tohle přidalo pár procent výkonu pro řazení intů, ale vedlo to k drastickému zpomalení v případě floatů. Výsledek nebyl nijak přesvědčivý a proto jsem tuhle optimalizaci nakonec nepoužil.
  4. To s přihládnutím k tomu, že radix sort musí udělat 5 průchodů celým polem znamená 2.2ns/7taktů CPU na element v každé iteraci. Z paměti je přitom třeba číst/zapisovat 4.8 GB/s (vysvětlení v tomto paperu, který hlásí, že jeden průchod jejich radix sortu, který dává 240M/s na čtyřech jádrech, trvá 13.2 taktu na element – pozor jejich čísla nesedí).
  5. V literatuře se uvádí, že radix sort nedělá sekvenční zápis v tom smyslu, že nezapisuje do jednoho místa, ale do mnoha, v tomto případě 256. To vede k neoptimálnímu využití propustnosti pamětí a může vést k neustálým výpadkům TLB. Moje experimenty na moderním hardwaru však nejsou přesvědčivé aspoň, co se týká TLB. Perf hlásí skoro nulu dTLB-load-misses a dTLB-store-misses, ale na druhou stranu cache miss se stále nějaké vyskytují (jako load-miss tak i store-miss).
  6. Se složitostí radix sortu je to o něco komplikovanější. Je pravda, že běží v čase O(kn), ale povaha onoho k je trochu matoucí. Když bych měl n různých čísel, nejmenší počet bitů potřebných k reprezentaci každého z nich je log<sub>2</sub>(n). Ono O(kn) je tedy vlastně jen O(n log(n)) v přestrojení.

Flattr this!

Posted in Algo, CPU, Paměť | 3 Comments

Dualismus hardwaru a softwaru, strojů a virtuálních strojů

René Descartes věřil v myšlenku dualismu, který se dá diletantsky vysvětlit tak, že tělo a duše jsou dvě různé vzájemně neslučitelné kategorie.

V podobném duchu se nese jiný dualismus, který je na rozdíl od Descarterových tvrzení méně metafyzický, ale zato je přítomný a hmatatelný. Oba pak mají jedno společné: Jde o falešnou dichotomii i když všichni v hloubi duši toužíme po opaku.

Pochopitelně mluvím o dualismu hardware a software, strojů a virtuálních strojů, kdy jeden dělá to, co to druhý nemůže, kdy „myslící“ program ovládá „hmotu“ stroje. Tato demarkace na první pohled působí správně a uspokojivě, jako kdyby to tak skutečně bylo. Když ale začnu problematiku zkoumat podrobněji a vezmu to oklikou přes historii CPU architektur, toto rozdělení, před chvílí ještě krystalicky jasné, najednou začíná mizet v mlhách.

Všechno, co vypadá jako neměnné pravdy vytesané do kamene, jsou jen rozhraní a úrovně abstrakce.

Už instrukční sada procesoru (ISA) není ten pomyslný pevný bod, na kterém stojí všechny ostatní želvy, ale také jen abstrakce. Téměř žádný dnešní procesor (možná s výjimkou některých jednoduchých RICSů nebo DSP) nevykonává instrukce tak, jak jsou za sebou vyskládány v binárce, ale dynamicky je v hardwaru překládá na jednodušší (RISCovější, i když tohle označení už moc neznamená) mikro operace (μops). Jedna instrukce se může rozpadnout na několik μops nebo naopak může být běžná sekvence několika instrukcí přeložena na jednu mikro operaci (tzv macro-op fusion)1. Situaci dále komplikuje mikrokód, kdy jsou některé instrukce procesorem přeloženy na malý program (například trigonometrické funkce jsou mikrokódovány, práce se denormálními floaty nebo GATHER instrukce na Haswellech, vylepšená v následujících generacích)14. Instrukce tedy není nutně nejmenší jednotka práce2 a často ani není přímo implementována v hardware.

Tato organizace není použita zbůhdarma, ale má dobrý důvod. Není totiž nutné plýtvat křemíkem na každou funkcionalitu, která je v programech používána jen okrajově. Místo toho je možné použít mikroarchitekturu s menším počet vysoce optimalizovaných funkčních jednotek, které jsou sdílené mezi instrukcemi3. Dále to zjednodušuje návrh pipeline, protože není třeba počítat s operacemi, které mají divoce rozdílné latence, a vede to k efektivnějším out-of-order jádrům, které můžou vykonávat různé μop (které trvají ±stejný počet taktů) různých CISC instrukcí najednou a využít tak dostupný ILP. Jako příklad můžu uvést stařičké VAXy. Ty například přímo podporovaly komplikovanou instrukci INDEX, která byla používána jen velice zřídka, a proto nikdo nestrávil čas optimalizací její hardwarové implementace. To došlo tak daleko, že když někdo napal stejnou funkcionalitu pomocí jednodušších operací, softwarový výsledek byl rychlejší než přímá implementace v křemíku. Víc o tom mluví Onur Mutlu v jedné ze svých skvělých přednášek o architektuře procesorů.

Jestliže instrukce jsou atomy, pak mikro operace jsou jejich kvarky.

Jako další příklad na spektru abstrakce může posloužit Transmeta. Ta v dobách své nedlouhé slávy produkovala čipy, které dokázaly rozběhat barokní x86 binárky. Sám hardware přitom neuměl nic z monstrózní x86 ISA, ale místo toho běhal na vlastní proprietární VLIW architektuře a x86 program byl za běhu překládán pomocí takzvaného Code Morphing software (CMS) do interní ISA4. Nepřekládal ale všechno. Zpočátku kód jen naivně interpretoval, když pak objevil opakovaně prováděný blok, rychle ho překompiloval z x86 do nativního VLIW formátu, když pak zjistil, že jde o skutečný hotspot, provedl důkladnou (a pomalou) kompilaci, jejímž výsledkem byl kvalitní a rychlý kód. Procesory od Transmety se tak v mnohém podobají JVM, jen s tím rozdílem, že jejich vstupem není Java bytekód, ale x86 binárka, cíl JITu není x86, ale interní nízkoúrovňový VLIW formát a hardware je specificky navržen pro tento účel.

VLIW instrukce se v mnohém podobají mikro operacím – jde o operace na velice nízké úrovni, které věrně opisují chování a mikroarchitekturu hardware5. Jsou prováděny ve statickém pořadí určeném kompilátorem (statically scheduled) a obnažují všechny interní detaily procesoru, jako jsou počty funkčních jednotek, délku pipeline a latence jednotlivých operací. To vede k problémům s přenositelností, kdy program zkompilovaný pro jeden model VLIW stroje nefunguje na novějším modelu, který má víc nebo rychlejší funkční jednotky.9 Když je však tento interní formát skryt pod jednou úrovní abstrakce (jako v případě x86 a dalších ne-VLIW ISA), může se mikroarchitektura libovolně měnit a vrstva o jednu úroveň výše se postará, aby byly všechny prostředky co nejlépe využity – nezáleží jestli jde o celkem naivní překlad v hardware, nebo dynamickou rekompilaci v software. Dokonce je teoreticky možné takto postavit stroj, který zvládá několik různých ISA najednou (třeba x86 a ARM).

V případě Transmety se interní formát měnil celkem znatelně – první generace procesorů pojmenovaná Crusoe byla široká 128 bitů a každá instrukce obsahovala 4 operace. Další generace byla dvakrát širší s dvakrát větším počtem operací v instrukci. Více informací je v oficiálním whitepaperu nebo ve dvou článcích o snahách reverzně vyinženýrovat mikroarchitekturu Transmety.

Velice podobně jako čipy už dlouho mrtvé Transmety pracují některé ARM čipy z rodiny Tegra od Nvidie založené na mikroarchitektuře Denver (nebo ruský Elbrus, když už jsem u toho). Na rozdíl od Transmety, která začala program v interpretru, Denver umí přímo vykonávat ARM kód na celkem základní úrovni. Hlavní síla čipů však spočívá v proprietární interní VLIW ISA, která je 7-wide (7 ops v instrukci). Program je zpočátku vykonáván v ARM módu s tím, že jakmile software objeví horký kus kódu, začne kompilovat. Rozpoznávání hotspotů nefunguje na úrovni základních bloků, ale na úrovni průchodů (trace), které jsou typicky tvořeny jednou iterací smyčky a můžou procházet mnoho základních bloků a dokonce i volání funkcí. Když je nalezena frekventovaná trace ARM instrukcí, je softwarově překompilována do VLIW ISA jako lineární sekvence operací, bez větvení a bez odboček. V mnohém tedy připomíná trasovací (tracing) JIT kompilátory, jako je třeba ten v PyPy (používaný i v jiných jazycích než Python), jen netrasuje nějakou formu bytekódu nebo operace interpretru (jako v případě meta-tracing JIT), ale přímo ARM instrukce.

Jako další příklad může sloužit společnost Azul, která dříve vyráběla čipy Vega specializované pro běh Java aplikací. Interně šlo o velice jednoduché a nepříliš rychlé in-order RISC procesory, které byly navrženy s ohledem na specifikaci Javy. Nepokoušely se však přímo vykonávat Javovský bytekód – lidé se o to dříve pokoušeli a zjistili, že je to velice špatný nápad. Na rozdíl od Teger a čipů Transmety měly Vegy vlastní operační systém, standardní knihovny a JVM zkompilované do nativního kódu, ale rozhraním pro uživatelské aplikace byl jedině Java bytekód. Šlo tedy o poměrně vysokou úroveň abstrakce, která návrhářům10 dala velkou míru volnosti jak společně navrhovat michroarchitekturu a JIT. Hardware se mohl mezi generacemi měnit (a také se měnil), aniž by to narušilo kompatibilitu pro uživatelské aplikace.

Úroveň abstrakce může být však i jinde, o čemž svědčí chystaný procesor Mill. Jde opět o VLIW staticaly scheduled, open pipeline stroj, který šířkou, kdy v jedné instrukci může být až 33 operací, překonává všechno ostatní. ISA je zamýšlena jako cíl pro kompilátor a bude přímo vykonána procesorem na nejnižší možné úrovni bez dalších vrstev abstrakce. Od ostatních VLIWů se Mill liší v tom, že není navrhován jako jedna architektura, ale rodina architektur s různými schopnostmi a hardwarovými prostředky.12

VLIW má problém v tom, že jde o příliš nízkou úroveň abstrakce – v podstatě představuje obnaženou mikroarchitekturu. Když se hardware změní, software se musí přizpůsobit. Mill se tohle snaží vyřešit tím, že programy nebudou distribuovány jako binárky přímo určené ke spuštění, ale budou na něco abstraktnější úrovni s tím, že se před spuštěním nebo během instalace staticky specializují pro konkrétní hardware. Na rozdíl od příkladů v minulých odstavcích by nemělo jít o nijak komplikovanou rekompilaci, ale celkem přímé mapování požadavků programu na prostředky hardware. Mill tak abstrahuje nejen možnosti daného stroje, ale i binární kódování instrukcí, které může být specifické pro každý model a přesto si zachovává výhody statického schedule, který nepotřebuje komplikovaný dekodér nebo out-of-order hardware.6

O Millu zatím nemůžu prohlašovat nic konkrétního, protože se ještě nedá koupit reálný křemík. Všechny informace pocházejí ze série velice detailních přednášek, které uspořádal Ivan Godard a ukazuje v nich všechny vlastnosti, kterými se Mill liší od běžných strojů a architektur.

Z předešlých odstavců lze vypozorovat jedno společné téma: všechno je jen abstrakce na určité úrovni. Od velice nízké jako v případě Millu, přes o něco vyšší (x86 a ARM ISA) až po relativně vysokou jako je Java bytekód. A to nemluvím o FPGA, které všechno tohle staví na hlavu, protože dovolí přímo programovat hardware. A víte jak Intel nedávno koupil Alteru? To znamená jediné: dočkáme se procesorů s integrovaný­mi FPGA.


Teď když jsem tu napsal přes 1000 slov o tom, že hranice mezi hardwarem a softwarem je mlhavější, než se na první pohled může zdát, nabízí se jedna zřejmá otázka: Neexistují nějaké principy psaní programů, které nějakým zásadním způsobem benefitují hardwarusoftwaru?

Zcela nepřekvapivě: Ano.

Jedna vlastnost, ze které profitují všichni, je předvídatelnost. Pokud se strojový kód chová předvídatelně, může hardware rozpoznat vzory jeho chování a pomoci. Sázka na předvídatelnost vedla k obohacení procesorů o branch prediction jednotky, které se snaží uhodnout zdali program skočí nebo propadne7, a branch target buffery, který předvídá cíle skoků na registr (typické při implementaci virtuálních metod), cache a prefetching.

Předvídatelnost pomáhá softwaru v tom, že dynamická prostředí, jako je JIT kompilátor v JVM, může specializovat kód pro běžné případy a ty, které nastanou jen zřídka nebo vůbec, ignorovat. Takto může devirtualizovat a inlinovat virtuální metody, což otevře dveře celé plejádě dalších optimalizací. Pro dynamické jazyky je předvídatelnost ještě důležitější, protože pokud má program rozumně statickou strukturu, virtuální stroj může objevit skryté třídy a typy argumentů funkcí a specializovat pro ně kód8. Předvídatelnost je také benefitem pro tracing JIT, které nahrávají (trasují) lineární průchod kódem, který je často prováděn, a čím méně odskoků a odboček v něm je, tím jednodušší má práci13. Stejně tak práce Maxime Chevalier-BoisvertBasic Block Versioning, ocení předvídatelnost tím, že bude potřebovat vygenerovat menší množství verzí základních bloků. Z předvídatelnosti také těží přístup založený na specializaci AST (jako je Truffle+Graal backend pro Jruby), protože (očividně) není třeba provádět takové množství specializací, což vede k menším a rychlejším několikaúrovňovým polymorfním inline cache.

Další dobrou vlastností pro všechny zúčastněné je malý a kompaktní kód. Na úrovni hardwaru je to proto, že se kód vejde do cache. Nejde jen o ty hlavní úrovně jako L3/L3/L1I, ale také L0, která na některých strojích obsahuje dekódované mikro operace. To pomůže v případech, kdy hlavním úzkým hrdlem je dekodér instrukcí. Pokud je však smyčka velice těsná (něco jako pár desítek μops na Haswell/Broad­well/Skylake x86) a celá se vejde do loop bufferu, hardware rozpozná, že si vystačí s tím, co je v loop bufferu a vůbec nebude komunikovat s cache a dekodérem. To může vést ke zrychlení programu a úsporám energie. Navíc nepředvídatelnému a nerozumně rozvětvenému kódu hrozí, že skočí na místo, které ještě není v cache a bude muset čekat.11

Virtuální stroje kompilující kód po metodách (jako typické JVM) preferují malé metody, protože jim to dává volnost v rozhodování, zdali danou metodu inlinovat nebo ne. Inlinování zvětšuje výsledný kód a proto, kdyby to JIT začal dělat příliš velkoryse, kód by nakynul a nemusel by se vejít do cache, což by vedlo ke znatelnému zpomalení. Inlinování je důležité, ale jen v rozumné míře. Když se to přežene, může to uškodit. Z toho důvodu JIT typicky vybírá kandidáty pro inlinování heuristicky a jedna z nich je velikost metody s tím, že velké metody mají jen malou šanci být inlinovány. Například JVM miluje malé metody a když se rozjede, začne je spekulativně inlinovat několik úrovní hluboko. Malé metody zároveň vedou k organizaci kódu s malými funkčními celky, které mají jasně definovanou roli, kdy každá funkce dělá jednu věc, což je obecně dobrá věc.


Nakonec to vypadá, že všechno jsou to jen virtuální stroje stojící na ramenou jiných virtuálních strojů. Dokonce i samotné Céčko, kterému se někdy říká přenositelný assembler definuje chování virtuálního stroje, který není nikde implementovaný v hardware, ale shodou náhod se na skutečný hardware dá docela pěkně napasovat. Všechny virtuální stroje se nakonec od sebe liší jen tloušťkou abstrakce, kterou je třeba překonat, aby program dokázal rozpohybovat reálné elektrony na reálném křemíku.


Dále k tématu

  1. Například poslední Intelí x86 kusy překládají sekvence některých instrukcí, jako cmp nebo aritmetické operace, následovaných skokem na jeden μop. Je tedy možné, že všechny instrukce, které řídí smyčku (inc, cmp, jXX) skončí jako jediný μop přímo implementovaný v hardware. Tohle je podle mě dobrý způsob jak utratit dostupný křemík, protože těsné smyčky jsou časté a program v nich stráví hodně času.
  2. I když některé jednoduché instrukce jsou přímo implementované v hardware. Jak je vidět z některých mikroarchitek­totických schémat, x86 CPU mají několik dekodérů z nichž většina dokáže dekódovat jen jednoduché instrukce, které jsou přeloženy na jeden μop a jeden, který překládá komplikovanější operace nebo pracuje s mikrokódem. Dekódování může pro komplikované CISC ISA představovat úzké hrdlo.
  3. Například aritmetické instrukce, které mají cíl nebe zdroj v paměti jsou přeloženy na jeden μop, který provede load/store a je neplánovaný na některou z dostupných load/store jednotek a další μop samotné aritmetické operace. Pro programátora viditelná architektura tak může být register-memory, ale mikroarchitektura je typu load-store.
  4. CMS naběhl ještě před startem operačního systému a byl to jediný kus softwaru napsaný přímo pro interní VLIW architekturu, všechno ostatní byl tradiční x86 od operačního systému až po uživatelské aplikace.
  5. Jenom s tím rozdílem, že VLIW přesně pasuje na dostupné funkční jednotky a celá dlouhá instrukce je vykonána paralelně, zatímco každý μop je na out-of-order strojích plánován a prováděn dynamicky. VLIW zachycuje statický instrukční paralelismus, μop na OOO dynamický.
  6. Mill se v mnohém podobá Itaniu, které částečně zapadá do rodiny VLIW strojů – každá instrukce je „dlouhá“ a obsahuje několik operací (v tomto případě se skupině operací říká bundle) a explicitní označení datových závislostí mezi operacemi. Pokud je mikroarchitektura dané implementace Itania dostatečně široká, stroj může vykonat pár bundlů paralelně a nemusí sám dynamicky zjišťovat závislosti mezi operacemi. Záměrem bylo dosáhnout většího ILP na novějších strojích bez nutnosti out-of-order hardwaru a se zachováním binární kompatibility.
  7. V současnosti branch predictor dosahuje úctyhodné přesnosti. Z paperu Prediction and the Performance of Interpreters je vidět, že procesor dokáže předpovídat skoky ve smyčce interpretru a svým způsobem pracuje na meta úrovni – nepředvídá skoky v kódu, ale v kódu, který interpretuje kód.
  8. Specializace se dá dotáhnout dál než skryté třídy/specializace kódu. Paper Adaptive Just-in-time Value Class Optimization se zabývá specializací kompozitních datových struktur.
  9. I některé ne-VLIW ISA odhalují interní detaily mikroarchitektury. Jde například o delay-slot.
  10. Mezi mě mimo jiné patřil Cliff Click, který měl o procesorech Vega skvělou přednášku, která kdysi, spolu s Performance Anxiety od Joshe Blocha, nasměrovala můj zájem směrem k hardwaru a procesorům.
  11. Ukazuje se, že některým serverovým úlohám nestačí instrukčí cache.
  12. Ale ani toto není nový nápad. Sám Ivan Godard říká, že Mill je rodina architektur ve stejném smyslu jako IBM System/360.
  13. Tady stojí za zmínku trace cache Pentia 4, která pracuje v duchu tracing JITu a ukládá lineární sekvenci dekódovaných mikro-operací.
  14. Podle knihy The Soul of a New Machine Tracy Kiddera stroje v sedmdesátých letech přecházely od přímé implementace instrukcí v hardware na mikrokód.

Flattr this!

Posted in CPU, VM | Leave a comment

Někdy je nejchytřejší nedělat nic chytrého (další kapitola nekonečného příběhu o optimalizaci)

Nedávno se na Twitteru objevil článek Fast Nearest Neighbor Queries in Haskell, který představuje novou knihovnu, která dokáže za pomoci velice chytrých algoritmů dělat nearest neighbor dotazy rychleji než jakýkoli jiný nástroj na planetě.

Nearest neighbor dotaz je takový, který má najít nejbližší bod z dané množiny bodů v určitém metrickém prostoru. Naivně se dá spočítat tak, že jeden po druhém projde všechny body s tím, že si pamatuje ten nejbližší a jeho vzdálenost1. Chytré řešení všechny body zaindexuje do stromové datové struktury, která nějakým rafinovaným způsobem rekurzivně rozdělí prostor do oblastí, a během hledání většinu těchto regionů přeskočí, protože je garantováno, že v nich obsažené body nemohou být blíže než současný kandidát. Zmíněná knihovna k tomuto účelu používá cover tree.

Článek rozhodně stojí za přečtení, už jen proto, že popisuje, jak psát rychlý numerický kód v Haskellu. Zcela nepřekvapivě jde hlavně o jednu věc: efektivní využití CPU cache ať už pomocí inlinování datových struktur, aby nebylo třeba nahánět zdivočelé pointery, nebo cache oblivious datových struktur jako třeba Van Emde Boas. Výsledky těchto snah jsou shrnuty v několika grafech s výsledky benchmarků.

A právě tady do příběhu vstupuji já.

Když jsem poprvé zíral na výsledná čísla, něco mi v nich nehrálo. Něco mi na nich nesedělo. Jeden benchmark pracoval se 100000 body o 384 dimenzích, ke každému hledal jednoho nejbližšího souseda měřeného euklidovskou vzdáleností a běžel 13.6 minuty. To by na mém stroji odpovídalo víc jak 2600 miliardám taktů procesoru. Naivní řešení by muselo projít všech 100000 bodů a každý porovnat se 100000 ostatními body a při každém porovnání pracovat s 384 prvky vektoru. Dohromady by bylo třeba 3840 miliard operací.

Skoro to vypadalo, že naivní přístup by nemusel být o moc pomalejší.

Bruteforcing for fun and profit

Sedl jsem si tedy ke stroji a za chvíli vyšvihl jednoduchý program v Céčku, který hledá nejbližší sousedy nejnaivnější možnou metodou.

double euclidean_distance_sq(double *vecs, int veclen, int a, int b) {
  double d = 0;
  for (int i = 0; i < veclen; i++) {
    double t = vecs[veclen * a + i] - vecs[veclen * b + i];
    d += t * t;
  }
  return d;
}

int nn(double *vecs, int veccnt, int veclen, int idx) {
  int nearest = -1;
  double dist = INFINITY;
  for (int i = 0; i < veccnt; i++) {
    double d = euclidean_distance_sq_vec_fl(vecs, veclen, idx, i);
    if (i != idx && d < dist) {
      nearest = i;
      dist = d;
    }
  }
  return nearest;
}

struct idxdist {
  int idx; // -1
  double dist; // INFINITY
};

int veccnt = atoi(argv[1]); // počet vektorů
int veclen = atoi(argv[2]); // délka jednoho vektoru

double *vecs = loadVectorsOrGenerateThemOrWhatever(veccnt, veclen);
struct idxdist *idxdists = initializeIdxsToMinusOneAndDistToMinusInfinity(veccnt);

for (int i = 0; i < veccnt; i++) {
  idxdists[i].idx = nn(vecs, veccnt, veclen, i);
}

Brnkačka. Program zkompilovaný pomocí gcc -O3 -funroll-loops běžel na 10× menších datech 36 vteřin. Aby byl srovnatelný s chytrým algoritmem, musel by běžet něco mezi 8 a 10 vteřinami.

Bylo na čase přitlačit na pilu, na vektorizovanou pilu. Nepodařilo se mi donutit GCC, aby automaticky vektorizoval smyčky a generoval SIMD instrukce7, a proto jsem se musel uchýlit k drastickým opatřením. Našel jsem seznam intrinsics – speciálních funkcí, které GCC přeloží na odpovídající instrukce9 – a ručně jsem vektorizoval funkci, která počítá euklidovskou vzdálenost. Nová verze provádí operace se čtyřmi doubly najednou2,8 a GCC tuto smyčku ještě odroluje, což přidá dalších pár procent výkonu6.

double euclidean_distance_sq_vec(double *vecs, int veclen, int a, int b) {
  __m256d d = _mm256_set1_pd(0.0);
  for (int i = 0; i < veclen; i+=4) {
    __m256d aa = _mm256_load_pd(vecs + veclen * a + i);
    __m256d bb = _mm256_load_pd(vecs + veclen * b + i);
    __m256d t = _mm256_sub_pd(aa, bb);
    d = _mm256_add_pd(d, _mm256_mul_pd(t, t));
  }

  __m256d sum = _mm256_hadd_pd(d, d);
  return ((double*)&sum)[0] + ((double*)&sum)[2];
}

Výsledek není na pohled nijak pěkný, ale mě nešlo o krásu, ale jen a pouze o surovou rychlost. A ta se zvýšila: Vektorizované verzi na výpočet stačí jen 20 vteřin. To je lepší, ale zdaleka ne ideální.

Jako další krok jsem zlepšil využití CPU cache. I když jsou vektory uložené v paměti jeden za druhým a procesor maskuje latenci pamětí načítáním dat dopředu (prefetch), limitujícím faktorem se stala propustnost paměti. Kód počítá vzdálenost jednoho vektoru se všemi ostatními. Onen jeden vektor si drží v cache, a všechny ostatní streamuje z DRAM. Každého se dotkne jednou a pak se hned přesune na další. Když začne porovnávat všechno s dalším vektorem, první vektory jsou už dávno vyhozené z cache a je třeba je opět streamovat z DRAM. To vede k mizernému využití cache.

Řešením je rozdělit pole vektorů do malých bloků, které se vejdou do cache, a počítat podobnost vždycky v jednom z těchto bloků. Když je velikost bloku 64, můžu porovnat 64*64 kombinací, ale z paměti stačí načíst jen 64+64 vektorů.3

int tile = 64;

for (int ti = 0; ti < veccnt; ti += tile) {
  for (int tj = 0; tj < veccnt; tj += tile) {
    int maxi = MIN(ti+tile, veccnt);
    for (int i = ti; i < maxi; i++) { // current vector index
      int maxj = MIN(tj+tile, veccnt);
      for (int j = tj; j < maxj; j++) {

        float d = euclidean_distance_sq_vec_fl(vecs, veclen, i, j);
        if (j != i && d < idxdists[i].dist) {
          idxdists[i].idx = j;
          idxdists[i].dist = d;
        }

      }
    }
  }
}

Kód obsahuje 4 vnořené smyčky, dvě vnější iterují po blocích (tile) a dvě vnitřní iterují obsahem bloků. Tohle uspořádání nevypadá na pohled příliš pěkně, ale plní svůj účel: program běží už jen 8.4 vteřiny a je tedy stejně rychlý jako chytrý algoritmus.

Ale to zdaleka nebylo všechno. V ten okamžik jsem se teprve dostal do ráže a v rukávu jsem měl připraveny ještě dva triky.

Dalším krokem bylo použití float místo double. Float potřebuje 4 bajty namísto osmi, z paměti tedy procesor dokáže protlačit 2× více dat, má 2× lepší využití cache a vektorové instrukce pracují s 8 floaty namísto 4 doublů.

float euclidean_distance_sq_vec_fl(float *vecs, int veclen, int a, int b) {
  __m256 d = _mm256_set1_ps(0.0);
  for (int i = 0; i < veclen; i+=8) {
    __m256 aa = _mm256_load_ps(vecs + veclen * a + i);
    __m256 bb = _mm256_load_ps(vecs + veclen * b + i);
    __m256 t = _mm256_sub_ps(aa, bb);
    d = _mm256_add_ps(d, _mm256_mul_ps(t, t));
  }

  d = _mm256_hadd_ps(d, d);
  d = _mm256_hadd_ps(d, d);
  return ((float*)&d)[0] + ((float*)&d)[4];
}

Tahle funkce je podobná té minulé, jen o něco málo ošklivější a 2× rychlejší. Program s těmito změnami běží jen 4.2 vteřiny – 2× rychleji než chytrý algoritmus.

Všechny tyto změny a úpravy přinesly skoro desetinásobnou akceleraci programu, který sám o sobě nebyl nijak strašně neefektivní. Ale samy o sobě jen optimalizovaly jedno-vláknový program a nesnažily se o nějakou razantní změnu. Poslední logickou metou byla paralelizace. Jak se ukázalo, tato změna byla nejjednodušší a nejrychlejší – stačilo před vnější smyčku přidat

#pragma omp parallel for

a zkompilovat s přepínačem -fopenmp4. Ze smyčky se rázem stala paralelní smyčka, která dokáže do posledního vytížit všechna jádra v systému. Program běžel 1.1 vteřinu. To je podle mě docela dobrý výsledek.

Pro zajímavost paralelní verze běží na mém netbooku s Atomem (2 jádra + hyperthreading, 1.66GHz) 17 vteřin, ±1000 vteřin na první generaci RaspberryPi (1 jádro, 700Mhz) a 51 vteřin na druhé generaci RasberryPi (4 jádra, 900MHz, NEON SIMD). Atom bohužel nepodporuje AVX a tak bylo třeba použít SSSE3, které pracuje jen se 128-bitovými vektory.

Shrnuto do jedné tabulky:

základní verze 36s
vektorizace 20s
lepší využití cache 8.4s
float namísto double 4.2s
paralelizace, OpenMP 1.1s

Kompletní program je k dispozici v tomto gistu.

Slovo na závěr

Na závěr musím dodat pár věcí, které zchladí všechny, kteří už-už vyrážejí slavit do ulic.

Autoři testovali onen chytrý algoritmus na Xeonu taktovaném na 2.8 GHz, takže jejich výsledky jsou o něco málo lepší, než se můžou zdát. Z prezentovaných čísel a grafů se dají těžko vyčíst nějaké absolutní hodnoty, ale vypadá to, že do testů zahrnují jak konstrukci cover tree, tak i jeho dotazování. Z jednoho grafu je vidět, že konstrukce nedominuje celkovému času, ale zabere něco jako 20%-30%. Z jiného se pak dá odvodit, že 50–90% všech přístupů do paměti skončilo jako cache-miss a kolem poloviny všech taktů procesor čekal na paměť a neměl co na práci, zatímco můj naivní program každý takt vykonal 2 instrukce a mnoho z nich pracovalo s osmi floaty najednou. Benchmarky jsem (věrný ideálům diletantství) nesnažil spustit na vlastním stroji, takže je možné, že mi mohlo uniknout něco zásadního.

Navíc: I kdyby byl onen cover tree v některých případech pomalejší, můžou existovat situace, kde bude kralovat5 – jako třeba menší počet dimenzí, jiná metrika nebo větší počet datových bodů. K tomu má stále výhodu parametričnosti, kdy můžu snadno změnit datové typy, reprezentaci a funkci vzdálenosti a nepotřebuje hledat několik bodů najednou, aby využil výhody rozdělení do bloků.

Závěr je následující:

Moderní procesory jsou neuvěřitelně rychlé a někdy hloupé řešení s dobrými konstantami může být rychlejší než chytrý algoritmus, který má na papíře lepší asymptoty. Neslavte však předčasně, vždycky nejdřív musíte vědět, čeho je možné dosáhnout.

Dále k tématu:


  1. Pokud mě zajímá n-nejbližších sousedů, pak je budu udržovat v haldě, která mi umožňuje rychle přistoupit k maximu.
  2. Vektor musí mít délku dělitelnou osmi a funkce neřeší, co se stane, když to neplatí. Není problém to dodělat přidáním malé smyčky nebo switche na konci. Dále je třeba mít na paměti, že _mm256_load_pd potřebuje, aby byla paměť zarovnaná na 32 bajtů, jinak segfaultne. Existuje také _mm256_loadu_pd, která načítá nezarovnanou paměť, ale je o něco pomalejší. V tomto případě je třeba místo malloc použít posix_memalign nebo podobnou funkci, která alokuje zarovnanou paměť.
  3. Velikost bloku musí být zvolena tak, aby se data vešla do nějaké úrovně cache. Pokud je tento parametr dimenzován pro L1 cache, je výsledek rychlejší než když je dimenzován pro L2 cache, ale může pojmout méně dat. Nastavení parametru pro sdílenou L3 cache je ještě pomalejší a má navíc ten problém, že L3 cache je většinou sdílená všemi jádry a každé z nich má k dispozici jen odpovídající část kapacity. Hodnota tohoto parametru je závislá na mnoha parametrech a proto ji není dobré nastavovat napevno.
  4. Více informací o OpenMP se dá najít v knize Programming on Parallel Machines
  5. Tady přichází ke slovu curse of dimensionality, která říká, že čím víc rozměrů má daný prostor, tím jsou všechny úhly bližší pravému úhlu a všechny vzdálenosti jsou si podobnější. To má za následek, že algoritmy jako cover tree a podobné fungují hůře a hůře, protože nemají jak prostor rozdělit na výrazně odlišné oblasti.
  6. I když to nemusí být tak úplně pravda na procesorech, které mají micro-op cache. Na ostatních strojích není nadměrná duplikace kódu také příliš vítaná. Když se program nevejde di instrukční cache, může to mít neblahé následky.
  7. Auto-vektorizace na první pohled vypadá jako triviální optimalizace, kterou kompilátor musí dát na první dobrou. Ale ve skutečnosti to ale není vůbec snadné a často je třeba kompilátoru pomoci. Viz: Auto-vectorization with gcc 4.7.
  8. Paralelismus vektorových instrukcí není to samé, co paralelismus out-of-order superskalárních CPU. OOO paralelismus se pokouší najít různé nezávislé instrukce, které může vykonat najednou, vektorové instrukce naproti tomu provedou stejnou operaci na několika datových položkách najednou. OOO je MIMD, vektorové instrukce jsou SIMD.
  9. Pro všechny, kdo nechtějí používat intrinsics, tu je vector extension GCC.

Flattr this!

Posted in Algo, CPU, Paměť | 2 Comments