0xDEADBEEF

RSS odkazy english edition
««« »»»

Bitonic sort

6. 7. 2021 #řazení #SIMD #D

Situace se má takhle: je potřeba seřadit obrovské množství malých polí fixní předem známé délky. Každé pole má pouze 20 intů, ale jsou jich jsou desítky milionů, možná stovky. Jak to udělat co nejrychleji? Neptejte se proč. Když situace nastane, je pozdě na otázky, záleží jen na rychlosti.

Rutina sort ze standardní knihovny jazyka D pro malá pole použije insertion sort2 a seřazení 20 intů trvá cca 210 ns. To bude moje referenční hodnota.

Jako první možnost se nabízí radixsort. Ten je neporazitelný pro velká pole intů nebo floatů (a s modifikací i pro data proměnlivé délky). Sestavování histogramů má ale určitou režii, která pro malá data zcela dominuje. Ve výsledku je radixsort mnohem pomalejší, než obyčejný sort. Tudy cesta nepovede.

Když dojdou rozumné varianty, nastane okamžik vyhrabávat kryptické znalosti doteď ležící ladem. Proč třeba nepoužít řadící sítě obecně a bitonic sort specificky? Četl jsem o nich už kdysi dávno a co si pamatuji, nepodávaly úplně nejhorší výsledky.

z wikipedie

Řadící síť je věc, která na vstupu dostane určitý fixní počet hodnot a na výstupu je vyprodukuje v seřazeném pořadí. Uvnitř je tvořena dráty vedoucí data směrem od vstupů k výstupům a komparátory, které porovnají dva dráty a prohodí jejich signály, pokud jsou v nesprávném pořadí. Jde jen o porovnání a cmov, žádný control flow a žádné podmíněné skoky. Pro seřazení n položek potřebují provést O(n (log n)^2) porovnání.

Šlo by jen o neškodnou akademickou zajímavost, kdyby se řadící sítě nedaly efektivně implementovat přímo v hardware nebo na SIMD architekturách. Grafické karty jsou k tomu ideální, jde o gigantické SIMD procesory, ale perfektně poslouží i SIMD na klasickém CPU.

Bitonic sorter je pak jedna konkrétní řadící síť, která se dá pohodlně implementovat na SIMD strojích.

Bitonic sorter z wikipedie

Červené krabice nakreslené vertikálně nad sebou můžou běžet paralelně jako série SIMD porovnání, permutací a maskování (shuffle, min, max, blend).

Čtyři prvky se seřadí tímto způsobem:

int4 sort4(int4 v) {
  // a <=> b, d <=> c
  auto v2 = shuffle(v, swapPairs);
  auto vmin = min4(v, v2);
  auto vmax = max4(v, v2);
  v = cast(int4) blend(cast(short8)vmin, cast(short8)vmax, 0b11000011);

  // a <=> c, b <=> d
  v2 = shuffle(v, swapHalves);
  vmin = min4(v, v2);
  vmax = max4(v, v2);
  v = cast(int4) blend(cast(short8)vmin, cast(short8)vmax, 0b11110000);

  // a <=> b, c <=> d
  v2 = shuffle(v, swapPairs);
  vmin = min4(v, v2);
  vmax = max4(v, v2);
  v = cast(int4) blend(cast(short8)vmin, cast(short8)vmax, 0b11001100);

  return v;
}

private enum swapHalves = 0b01001110;
private enum swapPairs  = 0b10110001;
private alias min4    = __builtin_ia32_pminud128;
private alias max4    = __builtin_ia32_pmaxud128;
private alias shuffle = __builtin_ia32_pshufd;
private alias blend   = __builtin_ia32_pblendw128;

Jde o přímou sekvence instrukcí, začne nahoře, skončí dole, žádné podmínky, žádné skoky, čas je vždy stejný, nezáleží na vstupních datech, výsledný vektorový registr obsahuje seřazené inty. Když se řadí dvojnásobek dat, počet kroků naroste o extra log(n), proto celková složitost n (log n)^2.

Asymptoticky jde o pomalejší algoritmus než obvyklé quicksorty, heapsorty a mergesorty, ale pro malé n je teoretický rozdíl zanedbatelný, prakticky to běží velice rychle a v případě, který mě zajímá, je bitonic sort mnohem rychlejší než vše ostatní.

Napsal jsem tedy bitonic sort specializovaný pro řazení 20 intů (zdroják tady). Používá jen SSE4 s 128-bit vektory, nikoli 256 bitové AVX2 z důvodu, že aktuální stroj (i5-2520M) nedisponuje AVX2 a AVX1 nemá velkou podporu pro integerové SIMD.

std.algorithm.sorting.sort210 ns10.50 ns/el
bitonic sort pro 20 intů41 ns2.05 ns/el
bitonic sort pro 16 intů17 ns1.06 ns/el

Velikost 20 je nejmenší blbé číslo. Pro 16 (jako mocninu dvou) stačí prostý bitonic sort - seřadit 4 čtveřice, pak sloučit dva páry čtveřic a nakonec sloučit jednu osmici - ale pro 4 inty navíc se musí dělat něco extra, jmenovitě vektorizovaný insertion sort. Výsledek není čistý bitonic sort, ale hybrid všeho, co se ukázalo jako nejrychlejší.

Drobné úpravy ve jménu rychlosti byly potřeba. Problémem se například ukázalo první řazení čtveřic (sort4 v ukázce nahoře). To má latenci kolem 8 taktů a je ho třeba provést 5x, což není málo. Naštěstí toho samého samého se dá dosáhnout i jinak a o něco rychleji. Když vezmu 4 vektory a paralelně je proženu bitonickou sítí nakreslenou v první obrázku nahoře. Výsledkem pak bude, že první elementy výsledných vektorů jsou vzájemně seřazené, to samé platí pro druhé, třetí a čtvrté. Pak stačí čtveřici transponovat a dostanu 4 seřazené vektory. To celé trvá pouhých 9 taktů. (Ve zdrojácích se o to stará funkce sort4x4 - min/max můžou běžet paralelně v jednom taktu, mají CPI propustnost 0.5, stejně jako intové punpck*).

Tyto 4 vektory dále řadím bitonic sortem jak je naznačeno na obrázku výše. Nakonec pak vektorizovaným insertion sortem přidám poslední 4 inty. Schéma použité sítě vypadá následovně.

Na začátku je vidět řazení čtyř čtveřic a následné transpozice, pak tělo bitonic sortu před nímž je nutné uvést některé vstupy do správného pořadí a nakonec diagonála insertion sortu.

Blbé je, že program natvrdo počítá jen s délku 20 a nic jiného neumí. V jazyce D by šlo napsat templatovanou funkci, která by specializovaný kód generovala na základě požadované délky v duchu ChipSortu. Výsledek by pak mohl být nejen velice svižný, ale i flexibilní a nabídnout to nejlepší z obou světů.


K tématu:


  1. Pro srovnání v Hoare’s Rebuttal and Bubble Sort’s Comeback hlásí s jejich variantou quicksortu časy 16.4 ns/element na poli 100k intů.
  2. Sledujte se mnou: sort volá quickSortImpl, ten volá shortSort, ten zavolá sort5 a pak provede insertion sort. sort5 je rutina bez smyček pro seřazení 5 položek.
píše k47 (@kaja47, k47)