Paralelní algoritmy
& architektury

Proč nás paměť vůbec brzdí — Memory Wall

Výkon procesorů rostl historicky mnohem rychleji než rychlost operační paměti RAM. CPU dnes umí v jednom taktu vykonat víc instrukcí, než mu paměť stihne dodat dat — tento rozevírající se nůžkový jev se nazývá memory wall (nebo von Neumann bottleneck). Rychlost paměti omezují dvě veličiny:

Statické vs. dynamické paměti

Adresování DRAM — řádek a sloupec

Adresa se dělí na RAS (row address strobe — index řádku) a CAS (column address strobe — index sloupce). Při čtení několika buněk ze stejného řádku stačí poslat jeden RAS a několik CAS → hardware sám zvýhodňuje sekvenční přístup v malých blocích.

Typy dynamických pamětí

• SDR SDRAM (Single Data Rate) — synchronní DRAM, taktovaná pomaleji než CPU.
• DDR SDRAM (Double Data Rate) — přenáší data na nástupné i sestupné hraně hodin → dvojnásobek dat za takt. Frekvence i takt zůstávají stejné, ale efektivně čteme ze dvou čipů.
• DDR2 → DDR5 — postupné zvyšování propustnosti širšími sběrnicemi a vyšším počtem paralelních čipů (banks). Latence se mění minimálně.

Cache — vyrovnávací paměti

Cache je malá rychlá SRAM mezi CPU a RAM. Velikost ~1/1000 RAM. Moderní procesory mají tři úrovně:

• L1 — nejmenší, nejrychlejší, na jádře. Typicky 32–64 KB.
• L2 — větší, sdílená v rámci jádra nebo dvojice jader. ~256 KB – 1 MB.
• L3 — sdílená napříč všemi jádry, řády MB.

Pravidlo: čím menší číslo, tím menší ale rychlejší cache.

Prefetching a cache line

CPU si proaktivně načítá data, o která ještě nepožádalo — odhadne z lokality program a načte celé bloky (cache line, typicky 64 bytů). Důsledek: čtení jednoho intu z RAM vždy přitáhne celou cache line s 16 sousedy zdarma.

Mapování cache — tři architektury

Write-through vs. write-back

• Write-through — zápis propisuje jak do cache, tak okamžitě do RAM. Bezpečné, ale pomalé.
• Write-back — zápis jen do cache; do RAM se zapíše až při vyřazení cache line. Rychlé, ale potřebuje protokol pro konzistenci (MESI).

⚠ Sekvenční vs. náhodný přístup — klíčový obrázek

Tohle profesor u zkoušky chce nakreslit:

• Sekvenční přístup: křivka taktů zůstává plochá a nízká — prefetcher má data připravená v cache.
• Náhodný přístup: jakmile dataset překročí kapacitu cache, počet taktů strmě stoupá — každý prvek vyvolá cache miss a CPU čeká stovky taktů na RAM.

Optimalizace z pohledu programátora

1. Sčítání matic po řádcích vs. po sloupcích

Matice jsou v C/C++ uloženy po řádcích. Iterace po řádcích = sekvenční přístup. Iterace po sloupcích = skok přes celý řádek za každým prvkem → totální cache miss. Rozdíl bývá řád.

// dobré: row-major iterace
for (int i = 0; i < n; i++)
    for (int j = 0; j < n; j++)
        C[i][j] = A[i][j] + B[i][j];

// špatné: sloupcová iterace skáče po n prvcích
for (int j = 0; j < n; j++)
    for (int i = 0; i < n; i++)
        C[i][j] = A[i][j] + B[i][j];

2. Násobení matic — transpozice nebo bloky

Klasické C = A·B potřebuje matici B procházet po sloupcích.

• Transpozice B: umožní obě matice číst sekvenčně, ale zvýší paměťové nároky.
• Blokové násobení (tiling): rozdělíme matice na malé bloky, které se vejdou do L1 cache, a každý blok výsledku spočítáme z bloků operandů. Nepotřebuje extra paměť a typicky vyhrává.

3. Critical stride — zákeřný hardwarový jev

4. Pole struktur vs. ukazatel na velká data

Když vyhledáváme podle klíče, nepotřebujeme tahat celé velké hodnoty do cache. Náhrada BigData data za BigData* data ve struktuře způsobí, že se do cache vejde víc klíčů a velká data se tahají až tehdy, když jsou opravdu potřeba.

CISC vs. RISC

• CISC (Complex Instruction Set) — bohatá instrukční sada, různé délky instrukcí, jedna instrukce může dělat víc věcí najednou. Klasické x86.
• RISC (Reduced Instruction Set) — málo, jednoduché, stejně dlouhé instrukce, lépe se zpracovávají v pipeline. Moderní x86 vnitřně CISC instrukce překládá na RISC mikroinstrukce.

Pipelining — zřetězené zpracování

Pipelining je forma paralelizace na úrovni CPU. Zpracování jedné instrukce se rozdělí na několik jednodušších kroků (jako výrobní linka):

Délka pipeline se liší architekturou: Pentium P5 = 5 fází, Core ~14, Netburst (extrémní) 20–31. Delší pipeline = vyšší frekvence, ale větší ztráta při špatném odhadu skoku.

Závislosti (hazardy) v pipeline

1. Datová závislost — instrukce čeká na výsledek předchozí, který ještě není hotový.
2. Závislost na zdrojích — dvě instrukce chtějí stejný registr / paměťové místo.
3. Závislost na podmínce (větvení) — průměrně každých 5–6 instrukcí, CPU neví, kterou větví jít.

Out-of-order & spekulativní zpracování

• Out-of-order execution — CPU (nebo už překladač) přehazuje pořadí instrukcí tak, aby se vyhnul blokacím datovými závislostmi.
• Register renaming — CPU flexibilně přejmenovává registry, aby uvolnil falešné závislosti.
• Spekulativní zpracování — když narazí na podmíněný skok, CPU uhádne větev a začne ji preventivně počítat. Pokud trefí, super; pokud ne, musí celou pipeline vyprázdnit (flush) a začít znovu — desítky taktů pryč.

Branch prediction — předvídání podmínek

CPU má speciální Branch Target Buffer (BTB) — malá cache (zhruba velká jako L1) ukládající historii skoků. Na základě stavového automatu (2bitový čítač) odhaduje pravděpodobnost, kam skok povede. Úspěšnost u dobře předvídatelného kódu ~99 %.

Superskalární zpracování

Pokud má program na sobě nezávislé instrukce, moderní CPU má víc pipeline (typicky 2–4) a dokáže provádět víc instrukcí současně. Aby to fungovalo, instrukce musí být dobře seřazené už od překladače.

Knihovna PAPI (Performance API) umožňuje vyčíst z procesoru statistiky o úspěšnosti predikce a využití pipeline.

Loop unrolling — rozbalení smyček

Každý cyklus má skrytou podmínku — kontrolu konce iterace. CPU musí v každém kole tipovat, jestli skok zpět nastane → úzké místo predikce.

// klasická smyčka — n podmínek
for (int i = 0; i < size; i++)
    sum += data[i];

// rozbalená 4×  — n/4 podmínek
for (int i = 0; i < size; i += 4) {
    sum += data[i];
    sum += data[i+1];
    sum += data[i+2];
    sum += data[i+3];
}

Provedení: ručně (přímo v kódu, případně C++ šablonami) nebo překladačem (gcc -funroll-loops). V testech ruční bývá lepší než překladačové.

Vektorizace (SIMD)

Procesory mají speciální širší registry (128 / 256 / 512 bitů) a nové sady instrukcí. Jedna instrukce (např. addps — add packed single) sečte najednou několik čísel. Architektuře se říká SIMD (Single Instruction, Multiple Data).

Kdo vektorizuje: typicky překladač s -mavx, -march=native apod. — programátor zasahuje jen pokud překladač zaváhá.

Co brání vektorizaci

1. Špatné zarovnání paměti — registry se nejlíp plní z adres zarovnaných na svou šířku. Nezarovnané čtení je dražší.
2. Aliasing — překladač neví, jestli se pointery v cyklu (např. b[i] += a[i]) v paměti nepřekrývají. Pak vektorizaci zakáže pro jistotu (řeší se kvalifikátorem __restrict).
3. Složité podmínky uvnitř cyklu — pokud se za jistých okolností opouští smyčka, nedá se snadno „rozbalit" do vektoru.

SMP / UMA — model sdílené paměti

Klasický vícejádrový procesor (notebook, server do ~64 jader). Všechny CPU jsou si rovny a mají stejně rychlý přístup do jedné globální RAM:

• UMA (Uniform Memory Access) — přístup k jakékoliv adrese trvá stejně.
• SMP (Symmetric Multiprocessor) — víc CPU, jedna RAM, jeden BUS.

Cache-coherence problem

Každé jádro má vlastní L1/L2 cache. Když dvě jádra načtou stejnou proměnnou X = 7 a jedno ji změní na 2, druhé o tom neví a počítá se zastaralou (stale) hodnotou. Hardware to musí řešit.

Existuje i alternativní update protocol, který místo zneplatnění rozesílá nové hodnoty — ale spotřebuje víc sběrnice. V praxi vyhrál invalidate.

⚠ False sharing

Obrana:

• Padding mezi proměnnými různých vláken — alokovat tak, aby každá proměnná byla na své cache line.
• Klauzule reduction v OpenMP — mezivýsledky v privátních proměnných každého vlákna, slučují se na konci.
• Mezivýsledky vláken neukládat blízko sebe.

U moderních CPU je dopad mírný (desítky procent), ne tak fatální jako dříve, ale stále se vyhneme.

Atomické instrukce

Místo složitých kritických sekcí lze pro jednoduché operace použít atomické instrukce. Jsou implementované přímo v procesoru a nelze je přerušit:

• CAS (compare-and-swap) — porovnej a vyměň.
• TS (test-and-set) — testuj a nastav.
• RMW (read-modify-write) — atomické přečtení, úprava, zápis.

Atomické instrukce mají nižší režii než kritická sekce, ale s rostoucím počtem vláken jejich cena roste (kontence na lince).

OpenMP — standard pro programování sdílené paměti

Bavíme se o vláknech (threads), ne procesech. OpenMP přidává pragmy nad C/C++/Fortran a překladač z toho udělá paralelní kód.

#pragma omp parallel for num_threads(4)
for (int i = 0; i < n; i++)
    a[i] = b[i] + c[i];

Klauzule pro proměnné

• shared — jedna sdílená kopie všech vláken (implicitní pro globální).
• private — každé vlákno má vlastní nezávislou kopii, neinicializovanou.
• firstprivate — private + počáteční hodnota se okopíruje z vnějšku.
• lastprivate — private + hodnota z poslední iterace se po skončení zapíše do původní proměnné.
• reduction(op:var) — privátní mezivýsledky se na konci sloučí asociativní operací. Eliminuje race condition i false sharing.

Scheduling — jak rozdělit iterace

Synchronizace

• #pragma omp critical — kritická sekce, jen jedno vlákno najednou.
• #pragma omp atomic — atomická aktualizace jedné proměnné (rychlejší než critical).
• #pragma omp barrier — bariéra: žádné vlákno nepokračuje, dokud nedorazí všechna.
• #pragma omp single / master — kód provede jen jedno / master vlákno.

// příklad: bezpečná suma přes redukci
double sum = 0.0;
#pragma omp parallel for reduction(+:sum)
for (int i = 0; i < n; i++)
    sum += a[i] * b[i];

Proč distribuovaná paměť

Sdílená paměť narazí na limit kolem ~64 procesorů — paměťový subsystém přestane stíhat všechny požadavky a cache coherence problém narůstá kvadraticky. Řešení: každý procesor má vlastní lokální RAM a procesory komunikují přes síť (interconnection network). Historicky: Cray T3D (1993).

Místo jednoho jádra běží nezávislé procesy (každý má vlastní adresový prostor). Komunikace probíhá výhradně message passing — předáváním zpráv. Standard: MPI.

Topologie sítí

Matematika komunikace

Doba odeslání zprávy se rozkládá na tři složky:

Packet routing s pipeliningem

Zpráva se posílá jako vláček paketů. Jakmile první paket projde routerem, ten ho hned posílá dál a přijímá další — paralelně. Vzdálenost mezi uzly (l) se tak stává téměř zanedbatelnou a celý vzorec se zjednoduší:

kde m je velikost zprávy. Tedy: startup penalizace + lineární přenos.

MPI — programovací model

MPI (Message Passing Interface) je knihovna a standard. Programátor explicitně volá Send/Recv. Místo vláken: nezávislé procesy (každý má svůj rank).

#include <mpi.h>

int main(int argc, char** argv) {
    MPI_Init(&argc, &argv);
    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);
    // ... práce podle ranku ...
    MPI_Finalize();
}

Blokující vs. neblokující komunikace

Blokující bez bufferu

MPI_Send pošle žádost a čeká, až druhý uzel oznámí připravenost přijímat. Pak teprve pošle data. Bezpečné, ale ⇒ riziko deadlocku.

Blokující s bufferem

Data se zkopírují do bufferu síťové karty a Send hned vrací. Žádné čekání → žádný deadlock, ale kopírování je extra režie a buffer se může přeplnit.

Neblokující

MPI_Isend, MPI_Irecv — síťové kartě se řekne „začni posílat", proces hned pokračuje (např. počítá nezávislé věci). Později MPI_Wait/MPI_Test zkontroluje dokončení. Nejefektivnější — překrývá komunikaci s výpočtem.

Komunikační operace — kreslit a popsat

Tyto chce profesor nakreslit a popsat:

1. Broadcast (one-to-all): jeden uzel rozešle identickou kopii dat všem ostatním.
2. Reduce (all-to-one): všichni pošlou data jednomu, cestou se asociativní operací (suma, min, max) slučují na síťových prvcích → strom redukce.
3. Scatter: jeden uzel rozkrájí velké pole a každému pošle jeho kus.
4. Gather: opak — všichni pošlou kousek jednomu, ten poskládá zpět.
5. All-to-all (personalized): každý posílá každému jinak specifická data — síťový ekvivalent transpozice matice. Nejnáročnější operace na síť.

Existují i kombinované varianty: MPI_Allreduce (reduce + broadcast výsledku všem), MPI_Allgather (gather + broadcast), MPI_Scan (paralelní prefix-sum).

Filozofie CPU vs. GPU

Důvod existence GPU: rendering 3D scén nemá datovou závislost mezi pixely/trojúhelníky → ideální pro masivní paralelizaci. Tato architektura se ukázala jako vynikající i pro maticové operace a neuronové sítě.

Hardwarový model — co umět namalovat

Vlákna, bloky, warpy, grid

• Vlákno (thread) — jednotka výpočtu. Vlákno = void funkce běžící na jednom CUDA core.
• Warp — skupina 32 vláken, která zpracovává vždy stejnou instrukci (SIMT — Single Instruction Multiple Threads). To je hardwarová jednotka.
• Blok (block) — softwarová skupina vláken (typicky 256, max 1024). Celý blok běží na jednom SM a vlákna v něm sdílí shared memory a mohou se synchronizovat.
• Grid — kolekce všech bloků jednoho kernelu. Bloky mezi sebou nemohou komunikovat ani se synchronizovat.

Divergence warpů

Co když warp narazí na if (cond) {A} else {B}, kde půlka vláken má cond=true?

1. Vlákna nemohou dělat různé věci najednou.
2. SM uspí polovinu warpu, druhá provede A.
3. Pak naopak — první polovina provede B, druhá čeká.
4. Doba běhu = doba A + doba B → ztráta výkonu na polovinu (nebo víc).

Proto se v CUDA kódu vyhýbáme rozsáhlým if/else uvnitř kernelu.

Paměťový model GPU

• Globální paměť — VRAM/HBM, viditelná všem vláknům, ~3 TB/s ale latence stovky taktů.
• Sdílená paměť (shared memory) — uvnitř SM, ~100× rychlejší než globální, viditelná jen v bloku. Programátor ji explicitně řídí.
• Registry — per-vlákno, nejrychlejší.
• Konstantní paměť — read-only, cachovaná.

Coalesced memory access — sloučené přístupy

Pokud 32 vláken warpu přistupuje na 32 sousedních adres globální paměti, GPU sloučí čtení do jediné transakce. Pokud přistupují roztroušeně (misaligned), vyžádá si to vícero transakcí a sběrnice se ucpe.

Memory banks ve sdílené paměti

Shared memory je rozdělena na banky. Pokud dvě vlákna sahají do stejné banky, vzniká bank conflict a přístupy se serializují. Programátor by měl zařídit, aby každé vlákno sahalo do jiné banky.

Programování v CUDA

__global__ void vectorAdd(float* a, float* b, float* c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) c[i] = a[i] + b[i];
}

// host (CPU):
cudaMalloc(&dA, n*sizeof(float));
cudaMemcpy(dA, hA, n*sizeof(float), cudaMemcpyHostToDevice);
// spuštění: 256 vláken na blok, ceil(n/256) bloků
vectorAdd<<<(n+255)/256, 256>>>(dA, dB, dC, n);
cudaMemcpy(hC, dC, n*sizeof(float), cudaMemcpyDeviceToHost);

CUDA podporuje C++ šablony, lambda funkce (částečně), printf z kernelu, ale ne všechny knihovny.

Pravidla pro efektivní kód na GPU

1. Minimalizovat přenos CPU ↔ GPU — PCIe je úzké hrdlo. Data nahraj na začátku, výsledek vrať na konci, mezi tím komunikuj minimálně.
2. Optimalizovat přístup do globální paměti — využívat coalesced access, kde to jde.
3. Omezit divergenci vláken — minimalizovat if/else uvnitř warpu, raději věci uspořádat tak, aby v jednom warpu byla stejná větev.
4. Zvolit vhodnou velikost bloku — typicky 128–512, mocniny dvou. Dotaz occupancy spočítá efektivitu.
5. Velké úlohy — GPU se vyplatí jen u objemných výpočtů. Pro malé úlohy režie přenosu spolyká celý zisk.

Pokud potřebujeme intenzivně komunikovat mezi CPU a GPU, existují hybridní karty se sdílenou pamětí (NVIDIA Grace Hopper).

Paralelní redukce — definice

Sekvenční implementace je triviální for cyklus a má složitost T_S(n) = Θ(n).

Stromová paralelní redukce

Paralelní verze využívá binární strom: v každé hladině polovina vláken sečte dvojici sousedů. Po log₂n krocích máme výsledek.

Složitost (n = p)

Nákladově optimální redukce

Klíč: méně procesorů (p < n). Postup:

1. Sekvenční fáze: každý z p procesorů zredukuje sekvenčně svůj blok n/p prvků → čas Θ(n/p).
2. Paralelní fáze: stromová redukce p mezivýsledků → čas Θ(log p).

V OpenMP přímou podporu má klauzule reduction(+:sum); v MPI MPI_Reduce a MPI_Allreduce.

Prefix sum — definice

Příklad: vstup [3, 7, 2, 5] → inkluzivní [3, 10, 12, 17], exkluzivní [0, 3, 10, 12].

Paralelní prefix sum

Vypadá jako sekvenční problém (každý prvek závisí na předchozím), ale jde paralelizovat stromově — stejný princip jako redukce, jen si uchováme i mezivýsledky. Složitost stejná jako redukce:

Algoritmus s n procesory není nákladově optimální — opět faktor log n nadbytek práce.

Nákladově optimální prefix sum

Podobně jako redukce: p < n, blok na proces.

1. Sekvenční fáze: každý procesor sekvenčně napočítá vnitřní prefix sum svého bloku A_k. Označme S_k hodnotu posledního prvku tohoto vnitřního prefixu.
2. Paralelní fáze: z posloupnosti S₁, …, S_p napočítáme exkluzivní prefix sum Σ_k — paralelně.
3. Korekce: každý procesor přičte ke svým vnitřním prefixům konstantu Σ_k.

Segmentovaný prefix sum

Trik: definujeme operaci ⊕ na párech (hodnota, příznak):

Operace je asociativní → můžeme přímo použít standardní paralelní prefix-sum algoritmus.

Kde se prefix-sum používá

• Filter / compaction — z masky (0/1) vytvořit kompaktní pole označených prvků.
• Radix sort — výpočet pozic přihrádek z histogramu (počty → offsety).
• CSR rowPointers — počty nenul v řádcích → ukazatele na začátky řádků.
• Alokace bufferů proměnné velikosti, dynamické fronty BFS.

Quicksort — sekvenční verze

Proč je paralelizace netriviální

• První volání = celé pole — máme jen 1 task. Až později se tasků zdvojnásobuje.
• Citlivost na pivot — pokud pivot rozdělí pole nerovnoměrně (např. 90:10), procesory dostanou nerovnoměrnou zátěž.
• Sekvenční quicksort je už dost rychlý — paralelní verze musí zaplatit režii pro převrácení.

Quicksort se sdílenou pamětí

Nesetříděná posloupnost je rozdělená mezi p vláken. Každé vlákno i drží podposloupnost.

1. Volba pivota — jedno vlákno (např. master) vybere pivota a uloží ho do sdílené proměnné. Všichni si ho přečtou.
2. Lokální partition — každé vlákno i rozdělí své prvky na S_i (menší) a L_i (větší).
3. Přeskládání do pomocného pole — pomocí dvou sdílených ukazatelů s = 0 a l = n−1:
   • smaller = atomicAdd(s, |Si|) — atomicky získá pozici pro své menší prvky (S roste zleva).
   • larger = atomicSubtract(l, |Li|) — atomicky získá pozici pro větší (L roste zprava).
4. Každé vlákno zkopíruje své S do smaller a své L do larger.
5. Pole se prohodí, vlákna si nově rozdělí prvky úměrně počtu a celý postup se opakuje.

Quicksort s distribuovanou pamětí

Místo zarážek S/L se použije prefix-sum:

1. Pivot — vyberou se na jednom procesu a broadcastem rozešle všem ostatním.
2. Lokální partition — každý proces i rozdělí své prvky na S_i a L_i.
3. Prefix-sum pro pozice v globálním poli: s_i = Σ_j<i |S_j|,   l_i = Σ_j<i |L_j|
4. Komunikace — proces i pošle své S na pozici s_i (sem může jít na jiný uzel) a L na pozici l_i. Typicky MPI_Alltoallv.
5. Posloupnost je teď ve dvou polovinách (S vlevo, L vpravo) — zopakuj rekurzivně pro každou polovinu na polovině procesů.

Hyperquicksort — řešení špatného pivota

Složitost paralelního quicksortu

Definice — radix sort

Princip jednoho kroku LSD — count + prefix-sum

1. Histogram (count): pro každou číslici j spočítám počet prvků s touto číslicí na aktuální pozici.
2. Prefix-sum z count dá offsety přihrádek: offset[j] = Σ_k<j count[k] To je pozice prvního prvku přihrádky j v pomocném poli.
3. Přeřazení (scatter): každý prvek se zařadí na pozici offset[jeho_číslice], offset se inkrementuje. Stejné číslice zachovají pořadí → stabilní třídění.
4. Opakuj pro každou pozici, od nejméně po nejvíce významnou.

Paralelizace: každý procesor napočítá lokální histogram. Kombinovaným prefix-sumem se z lokálních histogramů určí globální offsety. Pak každý procesor zapíše své prvky na správné pozice.

Radix sort se sdílenou pamětí — tři pokusy

Algoritmus 1 — naivní

Histogramy v poli count[j][pid] — index podle číslice, pak podle vlákna.

Algoritmus 2 — přerovnání indexů

Pole count[pid][j] — primární index je vlákno, takže různá vlákna píší do jiných cache lines.

• ✓ false sharing odstraněn
• ✗ špatná datová lokalita — v každé iteraci se prvky přehazují napříč celým polem, příští iteraci jiné vlákno přebírá data vzdálená v paměti.

Algoritmus 3 — MSD první krok + LSD rekurzivně

Klíčové vylepšení:

1. První krok provedu jako paralelní MSD — prvky se rozsejou do přihrádek podle nejvyšší číslice.
2. Přihrádky se rozdělí mezi procesory — žádný prvek se už dál nestěhuje mezi vlákny.
3. Každý procesor pak třídí jen svou přihrádku standardním LSD.

Naměřené výsledky (16M prvků, SGI Origin 2000)

20× zrychlení mezi Alg. 1 a Alg. 3 — z toho 80 % je odstranění false sharing, zbytek lepší cache lokalita.

Radix sort s distribuovanou pamětí

1. Každý uzel provede lokálně jednu iteraci radix sortu, v count[pid] má rozmezí přihrádek.
2. MPI_Allgather — všechny uzly získají kompletní pole count.
3. Z počtů se vypočítá bucketDistribution — celkové velikosti přihrádek a kdo dostane kterou (může se rozdělit i jedna přihrádka mezi dva uzly pro vyvážení).
4. MPI_Alltoallv — výměna dat: každý uzel pošle své prvky správným uzlům.

Třídící sítě — komparátory

Bitonic sort patří mezi třídící sítě — pevné struktury nezávislé na datech. Základní stavební prvek je komparátor:

Bitonická posloupnost

Bitonické rozdělení (bitonic split)

Klíčová věta o bitonickém rozdělení

Bitonic sort — síť komparátorů

Stavíme bitonické posloupnosti shora dolů a pak třídíme:

• Délka 1 → bitonická (triviálně).
• Dvě délky 2 (jedna ⊕ rostoucí, jedna ⊖ klesající) → bitonická délky 4.
• Dvě délky 4 → bitonická délky 8…
• Bitonická délky n → setříděná délky n pomocí bitonického merge (log n kroků).

Složitost

Proč je bitonic sort ideální pro GPU

• Pevná struktura — komparátory se mapují přesně, žádné dynamické rozhodování.
• Žádná divergence warpů — všechna vlákna provádějí stejnou instrukci (komparátor) na různých datech. Ideální SIMT.
• Implementace v sdílené paměti SM je přímočará (pole velikosti bloku).
• Nákladově není optimální (Θ(n log²n) celková práce vs. Θ(n log n) u quicksortu), ale konstanta je malá a hardware to vykompenzuje.

Blokové distribuované násobení matic

Pro matice A, B, C ∈ ℝ^n×n a p procesů (mocnina dvou) rozdělíme každou matici na √p × √p bloků. Každý blok C_i,j má spočítat jeden proces jako:

Přímočará verze a její problém

Pro p = 9 (tedy 3×3 bloků) vypadá první krok výpočtu:

Cannonův algoritmus — rotace bloků

Gaussova eliminační metoda (GEM)

Paralelní GEM — přímý chod

1. V kroku k vlákno k (pivotní řádek) vydělí svůj řádek pivotem a uloží do sdílené paměti.
2. Všechna ostatní vlákna i > k najednou a paralelně odečtou β_ik-násobek řádku k → jeden paralelní krok.
3. Na GPU: vynechat sdílenou paměť pivota a každé vlákno si dělí samo — odpadne čekání.

Problém zpětné substituce

Řešení — diagonalizace (Gauss-Jordan)

Místo trojúhelníkové matice nulujeme prvky nad i pod diagonálou. Výsledkem je diagonální matice → žádná zpětná substituce (x_i = d_i).

• Provedeme více aritmetických operací (řádově dvojnásobek).
• Ale veškerá práce je plně paralelizovatelná.
• Na superpočítači se vyplatí: čekání procesů je dražší než redundantní práce.

Pro řídké matice se klasický GEM nehodí (vyplnění — fill-in). Místo něj specializované řešiče (PARDISO, SuperLU) nebo iterativní metody (viz otázka 11–12).

Tridiagonální soustava

Soustava Ax = d, kde matice A má nenulové prvky jen na hlavní diagonále a dvou přilehlých:

kde c₁ = 0 a b_n = 0. Vznikají např. při diskretizaci jednodimenzionálních PDR (1D vedení tepla, struny).

Thomasův algoritmus — sekvenční

Cyklická redukce — paralelní trick

Odvození jednoho kroku

Bereme řádek i a jeho dva sousedy:

• c_i−1x_i−2 + a_i−1x_i−1 + b_i−1x_i = d_i−1
• c_ix_i−1 + a_ix_i + b_ix_i+1 = d_i
• c_i+1x_i + a_i+1x_i+1 + b_i+1x_i+2 = d_i+1

Zvolíme β_i = c_i/a_i−1 a γ_i = b_i/a_i+1. Od prostřední rovnice odečteme β_i×první a γ_i×třetí. Eliminujeme x_i−1 i x_i+1. Výsledek:

kde:

• c_i⁽¹⁾ = −β_i c_i−1
• a_i⁽¹⁾ = a_i − β_ib_i−1 − γ_ic_i+1
• b_i⁽¹⁾ = −γ_i b_i+1
• d_i⁽¹⁾ = d_i − β_id_i−1 − γ_id_i+1

Tuto úpravu lze provést na všech sudých řádcích nezávisle a paralelně. Vznikne tridiagonální soustava poloviční velikosti jen pro sudé neznámé x₂, x₄, …, x_n. Rekurzivně opakujeme.

Složitost cyklické redukce

Úplná cyklická redukce

Místo zpětné substituce (dalších log n kroků) provedeme redukci v obou směrech najednou:

1. V prvním kroku eliminujeme liché řádky → soustava jen pro sudé neznámé.
2. Současně eliminujeme sudé řádky → soustava jen pro liché neznámé.
3. Vzniknou dvě nezávislé tridiagonální soustavy poloviční velikosti.
4. Na každou aplikujeme stejný postup rekurzivně.
5. Po ⌊log₂ n⌋ krocích máme n/2 nezávislých 1- až 2-rovnicových soustav → vyřešíme triviálně paralelně.

Proč jsou grafy obtížné na paralelizaci

Grafy mají nepravidelnou strukturu (každý uzel má jiný počet sousedů) a typické algoritmy (BFS, DFS) dělají náhodné přístupy do paměti — cache misses na cache misses. Geniální trik celé přednášky:

Polookruh (semiring)

BFS — průchod do šířky

Algebraická formulace

Buď x vektor s 1 na pozici startovního uzlu a 0 jinde. Operace v booleovském polookruhu:

Výsledek y má 1 na pozicích uzlů, do nichž se v jednom kroku dá dostat. Opakovaným násobením A^T(∧.∨) šíříme frontu.

SSSP — nejkratší cesta z jednoho zdroje

Dijkstrův algoritmus — proč se nehodí

1. Inicializuj: d(start) = 0, ostatní ∞.
2. Vyber nenavštívený uzel s minimálním d(u).
3. Relaxuj všechny jeho sousedy.
4. Označ navštívený, opakuj.

Funguje jen pro nezáporné váhy. Složitost O((n+m) log n) s haldou.

Bellman-Ford v tropickém polookruhu

Bellman-Ford je iterativní — v každém kroku zkouší zkrátit cesty přes nějakou hranu. Funguje i pro záporné váhy (ale ne pro záporné cykly). Naivní složitost O(n · m).

Dijkstra vs. Bellman-Ford — srovnání

Shrnutí — co se paralelizuje a vyplatí