0xDEADBEEF

RSS odkazy english edition

Malé řazení ve vektoru #2 - phminpos

x86 nadělilo další dáreček. Čirou náhodou jsem se dozvěděl o instrukci phminposuw a čirou náhodou mi došlo, že se opět dá zneužít pro řazení. V tomto případě není tolik zajímavý finální program, jako cesta která k němu vedla a vše, co se na ní člověk naučil.

Zkratka phminposuw znamená packed horizontal minimum position unsigned word. Přehledně: Packed je opak scalar a znamená, že operace pracuje s celým vektorem a ne jen prvním prvkem. Horizontal znamená, že jde o horizontální redukci elementů jednoho vektoru, nikoli vertikální operaci mezi dvěma vektory. Unsigned word je důsledek historie x86, která začala jako 16-bit architektura a strojové slovo znamená 2 bajty.1 A konečně minimum position prozrazuje, co instrukce dělá: Vyprodukuje minimum vektoru a pozici tohoto minima. Něco jako následující funkce, jen v jedné instrukci implementované zcela v hardwaru.

ushort[8] phminposuw(ushort[8] src) {
  ushort index = 0;
  ushort min = src[0];

  foreach (i; 1 .. 8) {
    if (src[i] < min) {
      index = i;
      min   = src[i];
    }
  }

  return [min, index, 0, 0, 0, 0, 0, 0];
}

Nejjednodušší způsob, jak ji použít pro seřazení vektoru, představuje smyčka, která najde minimum, odloží ho stranou, vymaskuje jeho pozici za maximální hodnotu 0xffff a proces opakuje celkem osmkrát.

void minposSort(ref ushort[8] vec) {
  ushort8 v = unalignedLoad(vec);
  ushort[8] res = void;

  static foreach (i; 0 .. 8) {{
    ushort8 minAndPos = hminposuw(v);
    res[i]     = minAndPos.array[0];
    ushort pos = minAndPos.array[1] & 0b111;
    v.array[pos] = cast(ushort) 0xffff;
  }}

  vec = res;
}

alias unalignedLoad = (vec) => cast(ushort8) __builtin_ia32_loaddqu(cast(char*) vec.ptr);
alias hminposuw     = (vec) => cast(ushort8) __builtin_ia32_phminposuw128(cast(short8) vec);

Za pozornost stojí maskování 0b111. Jde o to, že kompilátor (GDC 10.2) neví nebo neumí využít znalost, že pos se vždy nachází v rozsahu 0-7, což představuje validní index do vektoru short8. Bez téhle masky musí do každé smyčky přidat kontrolu, jestli index nepřestřelil, což stojí čas a o čas tu jde především. Přesto je výsledek pomalý. Mnohokrát pomalejší než by měl být.

Problém spočívá v tom, že řádek

v.array[pos] = cast(ushort) 0xffff;

nemá žádný přímý protějšek v SIMD ISA. (Kompiloval jsem pro nepříliš nový laptop, který je vyzbrojen jen SSE4 a AVX.) Když kompilátor neví, jak to efektivně přeložit, vypomáhá si pamětí. Na začátku každé iterace načte vektor z paměti, na konci ho uloží zpátky a přepíše index pos a i když zásah do L1 cache na použitém HW trvá 4 takty, jde o 4 takty v řetězu závislostí. OOO nemůže spekulovat kolem a musí čekat, až se data vrátí z L1$.

Problém se SIMD se dá vyřešit jen nsazením více SIMD, jak jinak.

  // ...

  ushort8 mask = [0, 1, 2, 3, 4, 5, 6, 7];
  static foreach (i; 0 .. 8) {{

    // ...

    //v.array[pos] = cast(ushort) 0xffff;
    v |= onesMask(mask, pos);
  }

alias onesMask = (mask, pos) {
  ushort8 v = pos; // broadcast
  return cast(ushort8) __builtin_ia32_pcmpeqw128(cast(short8) mask, cast(short8) v);
};

Ve funkci onesMask se pozice broadcastuje do vektoru, porovná se s poziční maskou mask, výsledek je vektor nul, kde jen ushort se správným indexem je plný jedniček (0xffff). Pak se v jen zoruje a je to. Příslušná pozice vymaskována.

Tohle je konečně dostatečně rychlé, aby se to vyplatilo měřit. Přesto v porovnání s verzemi z minulého článku nic moc.

(Sandy Bridge i5-2520M)
doNothing          2.29 ns/vector
minposSort        20.44 ns/vector
---
pcmpstrSort       13.39 ns/vector
vecPrefixSort     11.67 ns/vector

doNothing je funkce, která nic nedělá a měří jen režii volání funkce. Přesvědčil jsem se, že měří, co si myslím, že měří - výsledná binárka obsahuje příslušnou call instrukci, nedojde k inlinování.

void doNothing(ref ushort[8] arr) {
  pragma(inline, false);
  asm { ""; }
}

Ještě ale nejsem u konce sil. Dá se to celé napsat jinak a místo maskování odečítat.

void minposSortSub(ref ushort[8] vec) {
  ushort8 v = unalignedLoad(vec);
  ushort[8] res = void;
  ushort sub = 0;

  static foreach (i; 0 .. 8) {{
    ushort8 minAndPos = hminposuw(v);
    ushort minval = minAndPos.array[0];

    res[i] = cast(ushort) (minval+sub);
    v     -= cast(ushort) (minval+1);
    sub   += cast(ushort) (minval+1);
  }}

  vec = res;
}

Tato verze najde minimum, pozice se zahodí, minimum+1 se broadcastuje do vektoru, ten se odečte od toho vstupního (v) a tím pádem všechny hodnoty poklesnou. Současné minimum přeteče do maxima, druhý nejmenší element se stane novým minimum a smyčka může pokračovat. Vypadá po poněkud krkolomně, nefunguje pro duplicity, ale je to rychlejší. Pravda, jen marginálně, ale tady se hraje o jednotlivé takty.

minposSort        20.44 ns/vector
minposSortSub     19.98 ns/vector

Půl nanosekundy znamená asi jeden a půl taktu na mém stroji.

Ok, pohled do binárky odhalí následující řetěz vzájemně závislých instrukcí.

vphminposuw // (5) hminposuw(v)
vpextrw     // (2) minAndPos.array[0] (vector -> scalar)
add         // (1) minval+1
vmovd       // (1) (scalar -> vector)
vpshufb     // (1) broadcast
vpsubw      // (1) v -= ...

Podle Agnera tato sekvence na Sandy Bridge CPU trvá přinejlepším 11 taktů. Za pozornost stojí, že přechází z vektorových registrů do těch normálních, aby tam se tam přičetla jednička a pak zase putuje zpátky. To není nutné, SIMD taky umí sčítat.

void minposSortSubVec(ref ushort[8] vec) {
  ushort8 v = unalignedLoad(vec);
  ushort[8] res = void;
  ushort sub = 0;
  ubyte16 mask = [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1];

  foreach (i; 0 .. 8) {{
    ushort8 minAndPos = hminposuw(v);

    auto mins = shuffle(minAndPos, mask);
    v  -= (mins + 1);

    ushort minval = minAndPos.array[0];
    res[i] = cast(ushort) (minval+sub);
    sub   += cast(ushort) (minval+1);
  }}

  vec = res;
}

alias shuffle = (vec, mask) => cast(ushort8) __builtin_ia32_pshufb128(cast(ubyte16) vec, mask);

Program na první pohled vypadá celkem podobně, přesto ve výsledné binárce se udála velká změna. Řetěz závislých instrukcí ve smyčce se značně smrskl na pouhých 8 taktů.

vphminposuw // (5) hminposuw(v)
vpshufb     // (1) shuffle(minAndPos, mask)
vpxor       // (1) x = -mins-1
vpaddw      // (1) v += x

Třetí řádek xoruje se samými jedničkami a jde o způsob, jak provést a ^ 0xffff = -a-1 bez potřeby immediate.

Extrakce a zápis výsledku pořád musí dělat vpextrw, ale na rozdíl od předchozí varianty se nenachází v řetězci závislých operací. Nezávisí na ničem dalším a nic nezávisí na něm a může běžet paralelně bokem. Hardware má dost paralelních zdrojů, aby se o to postaral. Čas je v tomto případě limitovaný délkou řetězce závislých instrukcí.

minposSort        20.44 ns/vector
minposSortSub     19.98 ns/vector
minposSortSubVec  15.74 ns/vector

Několikrát to tu už zaznělo: Jde o plně sekvenční algoritmus, kdy jedna iterace smyčky plně závisí na té předchozí. Hardware nemá možnost dělat paralelně nic kromě malých vedlejších větví.

Tady může pomoct trik, o němž jsem poprvé četl tady v kontextu zrychlení merge sortu. Ten je také plně datově závislý a iterace nemůže pokročit dokud neskončila ta předchozí.

Trik spočívá v tom, že pole má dva konce a můžu ho zpracovávat stejným sekvenčním programem z obou konců najednou. Vzniknou tak dva proudy instrukcí, které jsou na sobě vzájemně nezávislé a hardware dostane příležitost je provést paralelně.2

void minposSort2(ref ushort[8] vec) {
  pragma(inline, false);

  ushort8 pos = unalignedLoad(vec);
  ushort8 neg = ~pos;
  ushort[8] res = void;
  ushort8 mask = [0,1,2,3,4,5,6,7];

  static foreach (i; 0 .. 4) {{
    ushort8 minAndPos = hminposuw(pos);
    ushort8 maxAndPos = hminposuw(neg);

    res[i]        = minAndPos.array[0];
    ushort minpos = minAndPos.array[1] & 0b111;
    pos |= onesMask(mask, minpos);

    res[7-i]      = ~maxAndPos.array[0];
    ushort maxpos = maxAndPos.array[1] & 0b111;
    neg |= onesMask(mask, maxpos);
  }}

  vec = res;
}

Všechny varianty vypadají podobně. Vstupní vektor je v proměnné pos, jeho negace pro iteraci pozpátku v neg. Počet iterací se snížil na polovinu, tělo smyčky obsahuje krok oběma směry, jeden funguje jak je nastíněno výše, druhý musí výsledek negovat. Maximum je stejné jako minimum negované proměnné.

                  jeden konec        oba konce
minposSort        20.44 ns/vector    15.79 ns/vector
minposSortSub     19.98 ns/vector    13.89 ns/vector
minposSortSubVec  15.74 ns/vector    11.61 ns/vector

Výsledek je skutečně rychlejší. Hardware, i když jde o starý model CPU, má rezervy v paralelismu a takhle je může dát k dobru. Na novějších mikroarchitekturách s větší mírou instrukčního paralelismu, propast mezi sekvenčním programem a dostupnými prostředky bude představovat větší problém. Tohle není jako paralelizace s pomocí vláken. Když program běží jen v jednom vlákně, na ostatních procesorech může běžet něco jiného. Tady jde o zdroj, který když nebude využitý, zůstane ležet ladem.

Ok, to všechno je zajímavé. Na datových závislostech záleží, hlavně na out-of-order architekturách, které dnes pohání všechno od sledovacích zařízeních ve vaší kapse po datová centra.

Ale jak probíhá měření? Měří, co si myslím, že má měřit?

static foreach (f; AliasSeq!(doNothing, minposSortSubVec, minposSortSubVec2, pcmpstrSort, vecPrefixSort)) {{
  auto timer = StopWatch(AutoStart.yes);
  foreach (i; 1 .. iters) {
    arr = src; // zkopíruje se, co se bude řadit
    f(arr);    // seřadí se, ale na výsledku nic nezávisí
  }
  writeln(fullyQualifiedName!f, " ", timer.peek.total!"nsecs" / double(iters), " ns/vector");
}}

Program skutečně volá, co volat má, díval jsem se do binárky, problém je v řádku arr = src. Zkopíruje testovací vektor, pak ho seřadí, ale výsledek se zahodí a nic na něm pak dál nezávisí. Dostatečně agresivní hardware může spekulovat dost napřed nato, aby paralelně prováděl několik řadících funkcí paralelně. Každá může mít dlouhý řetězec závislostí, ale procesor jich pár překryje.

Když zakomentuji řádek arr = src, najednou do testovacího programu vnesu datovou závislost. Jedno řazení nemůže začít dokud to předchozí skončilo. Obrázek pak vypadá následovně.

                    nezávislé          závislé
doNothing            2.26 ns/vector     1.95 ns/vector

minposSort          20.44 ns/vector    39.27 ns/vector
minposSortSub       19.98 ns/vector    40.52 ns/vector
minposSortSubVec    15.74 ns/vector    36.08 ns/vector

minposSort2         15.79 ns/vector    27.23 ns/vector
minposSortSub2      13.89 ns/vector    26.00 ns/vector
minposSortSubVec2   11.61 ns/vector    22.09 ns/vector

pcmpstrSort         13.39 ns/vector    23.27 ns/vector
vecPrefixSort       11.67 ns/vector    16.47 ns/vector

bitonicSort          5.75 ns/vector     9.95 ns/vector

Při nezávislém testování, řazení s hminposuw zrychlí asi dvakrát. vecPrefixSort z minulého článku zrychlí méně, přesto je nejrychlejší. Důvod je prostý: krátké sekvence závislých instrukcí a velký potenciál pro interní ILP. I když běží jen jedna funkce v jeden okamžik, CPU provede hodně práce paralelně.

Ani jeden přístup není správný nebo špatný. Záleží jen, co chci měřit, latenci nebo propustnost. V různých situacích záleží na různých pohledech, jestli měřená funkce představuje úzké sekvenční hrdlo nebo zdali jde o propustnost.

Za pozornost ještě stojí poslední řádek: bitonic sort je rychlejší než všechny ostatní varianty. Jde o sekvenci 31 jednoduchých SIMD instrukcí zcela bez skoků. Každá instukce je závislá na té předchozí, délka řetězu závislostí je 31 taktů, což na 3.2 GHz stroji dává těch ~10 ns.


  1. Zajímavé je, že artefakt historie žije dál v GPU architekturách intelu i AMD. 32 bitů se na obou označuje jako dword.
  2. Oboustrané zpracování je zmíněné taky v tomhle excelentním článku
píše k47 (@kaja47, k47)