Motivace
PA081: Programování numerických výpočtů
8. Lineární algebra rychle
Vyrovnávací paměti
Blokové algoritmy
Aleš Křenek
LU
dekompozice
Vektorové instrukce
jaro 2012
Motivace
► proč má smysl optimalizovat právě algoritmy lineární algebry?
► frekventované problémy
► řešení lineárních rovnic, včetně velkých
► řešení nelineárních rovnic
► optimalizace
► metoda konečných prvků
Motivace
proč má smysl optimalizovat právě algoritmy lineární algebry?
frekventované problémy
► řešení lineárních rovnic, včetně velkých
► řešení nelineárních rovnic
► optimalizace
► metoda konečných prvků
cíle této přednášky
► zdánlivě jednoduché algoritmy (násobení matic) lze výrazně vylepšit
► problematika se stále vyvíjí
► naznačíme směry, kterými se lze vydat
► rafinované algoritmy už jednou někdo vymyslel a implementoval
► nemá smysl učit se je nazpaměť
► stačí o nich vědět a rozumět základní myšlence
4 □ M0M1 M 1 » 1
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
■O0.O
2/18
Vyrovnávací paměti
hlavní paměť je řádově pomalejší než procesor
► rychlou paměť v plné velikosti je technicky/finančně nemožné vyrobit
hierarchie vyrovnávacích pamětí (cache) pro Intel Nehalem/Westmere
► LI - plná rychlost, 32 kB instrukce, 32 kB data/jádro
► L2 - latence 10 cyklů, 256kB/jádro
► L3 - latence 40 cyklů, 8 MB/procesor (sdílená mezi jádry)
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Vyrovnávací paměti
Organizace přístupu
► řádky (cache-line)
► typicky 64B (16x float, 8x double)
► přímo mapovaná cache
► blok dat z hlavní paměti má dán jeden řádek v cache
(adresa/velikost řádku) % počet řádků
► snadno dojde ke kolizím, např. pro 32 kB
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
ouble pole[100][512][8];
(i=0; i<100; a += pol e [i] [0] [0];
řešením je deklarace pol e [100] [513] [8]
■O0.O
Vyrovnávací paměti
Organizace přístupu
plně asociativní
► blok dat z hlavní paměti může být umístěn v kterémkoli řádku
► implementačně náročné, ve standardních CPU se nepoužívá n-cestná asociativní
► kompromisní řešení
► sady po n řádcích
► blok z hlavní paměti může skončit v kterémkoli řádku sady
► Intel Nehalem/Westmere: n = 8
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Vyrovnávací paměti
Praktické zásady
optimalizovat pro všechny úrovně využít celý řádek
► např. může se vyplatit transponovat matici
dovolit procesoru/kompilátoru současně přenášet data a počítat
► dobře čitelný kód bez možných vedlejších efektů zabránit předčasným kolizím v jedné sadě řádků měřit dosažený výkon (flop/s)
atd. ...
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Vyrovnávací paměti
Praktické zásady
optimalizovat pro všechny úrovně využít celý řádek
► např. může se vyplatit transponovat matici
dovolit procesoru/kompilátoru současně přenášet data a počítat
► dobře čitelný kód bez možných vedlejších efektů zabránit předčasným kolizím v jedné sadě řádků měřit dosažený výkon (flop/s)
atd. ...
aplikovat jen pro skutečně kritické sekce kódu používat speciální optimalizované knihovny, kde to jde
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
6/18
Blokové algoritmy
► zjednodušený model jen s Ll
► násobení vektoru matici y := y + Ax
načti x do cache načti y do cache for i = 1.....n
načti řádek A;* do cache
f o r j = 1.....n
y t =yi + AíjXj
zapiš y do hlavní paměti
► počet přístupů do pomalé paměti m = 3n + n2 *■ počet aritmetických operací / = 2n2
*■ efektivita f Im, « 2
► algoritmus je limitovaný rychlostí paměti
► pomůže recyklace řádků cache - vyplatí se ukládat spojitě, ne jako sloupce matice (v C)
► jinak se s tím nedá nic moc dělat
Blokové algoritmy
Maticové násobení
► násobení matic C = C + AB
for i = 1.....n
načti řádek A;* do cache for i = 1.....n
načti Cíj do cache
načti sloupec B for k = 1.....r,
do cache
Cíj = Cíj + AikBkj
zapiš Cíj do hlavní paměti
přístupy do paměti m = n2 4 aritmetické operace / = 2n3 efektivita opět « 2
2nl
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
8/18
Blokové algoritmy
Maticové násobení
► matice rozdělíme na N2 bloků velikosti b xb,b = n/N
do cache
for i = 1.....JV
for j = l.....N
načti blok C for k = 1,..".,N
načti blok Ajj; do cache načti blok By do cache
zapiš blok C j.- do cache
přístupy do paměti
► čtení bloků A: N3 (n/N)2
► dtto B: N * n2
► čtení a zápis C: 2n2
Cjj + AjfcBfc; (maticové násobení)
■ N * n2
efektivita f /m
Ire (2JV+2)n2
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
9/18
Blokové algoritmy
Maticové násobení
PA081: Programování numerických
výpočtů
A. Křenek
Motivace
velikost reálne LI cache je 32 kB,
Vyrovnávací paměti
► do LI se vejdou 3 matice velikosti 36 x 36 (v doubí e)
► žádoucí aplikovat rekurzivně i pro L2 a L3
► implementováno v optimalizovaných knihovnách
► násobení matic je na současných procesorech jeden z nej efektivnějších algoritmů
► redukce problému na mat. násobení při řádovém zachování počtu operací téměř vždy vede k výraznému zrychlení
L2 cache je lOx pomalejší procesor je dobře využit
► i při využití vektorových instrukcí
► proto je cache právě tak velká
Blokové algoritmy
BIAS
Asymptotická složitost
LU
dekompozice
■O0.O
Knihovna BLAS
Basic Linear Algebra Subprograms, http://www.netli b.org/blas základní operace, např.
► DOT:xy
► AXPY:y + Ax
► NRM2: |x|2
► TRMV: Ax
► TRSV:A_1x
► G EMM: fiC + oíAB
verze pro typy f 1 oat,doubí e,compl ex
specializované varianty pro symetrické, trojúhelníkové, a některé řídké matice
de-facto standardizované rozhraní
implementace vyladěné pro konkrétní typy CPU
využívané v knihovnách vyšší úrovně (např. LU dekompozice)
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
11/18
Asymptotická složitost
► časová složitost násobení matic je 0{n3)
PA081: Programování numerických
výpočtů
A. Křenek
Vyrovnávací
paměti
Blokové algoritmy
Asymptotická složitost
LU
dekompozice
Vektorové instrukce
Asymptotická složitost
► časová složitost násobení matic je 0(n3)
► Strassen (1969): 0(n281)
► rozdělení matic na 2 x 2 bloky, potom
Mi = (Au +A22)(Bii +B22)
M2 = (A21 + A22)Bii
M3 = An(Bi2 +B22)
M4 = A22(B2i - B11)
M5 = (Au + Ai2)B22
M6 = (A2i -Aii)(Bn +B12)
M7 = (A12 -A22MB21 +B22)
C11 = Mi + C12 = M3 + C2i = M2 + C22 = Mi -
M4 - M5 + M7
M5
M4
M2 + M3 + M6
rekurzivní aplikace, počet operací T(n) T{n) = 7T(n/2) + 18(n/2)2 = 0(nlo&7)
0(n
2.811
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
□ 4 :
12/18
Asymptotická složitost
► časová složitost násobení matic je 0(n3)
► Strassen (1969): 0(n281)
► rozdělení matic na 2 x 2 bloky, potom
Mi = (Au +A22)(Bii +B22)
M2 = (A21 + A22)Bii
M3 = An(Bi2 +B22)
M4 = A22(B2i - B11)
M5 = (Au + Ai2)B22
M6 = (A2i -Au)(Bn +B12)
M7 = (A12 -A22)(B2i +B22)
C11 = Mi + M4 - M5 + M7
C12 = M3 + M5
C2i = M2 + M4
C22 = Mi - M2 + M3 + M6
rekurzivní aplikace, počet operací T(n)
T{n) = 7T(n/2) + 18(n/2)2 = 0(nlo&7) Coppersmith-Winograd (1987): 0(n2376)
0(n
2.811
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
12/18
LU dekompozice
► rekurzivní formulace algoritmu
/ A,
skalární a bloková verze
a2i
»12
= a2i/ofn zůstává
A22 = A22 - a2iai2
podstatná část operací je násobení matic takto implementuje knihovna LAPACK
□     s   4 :
Vektorové instrukce
13/18
LU dekompozice
► rekurzivní formulace algoritmu
/ A,
skalární a bloková verze
a2i
»12
= a2i/ofn zůstává
A22 = A22 - a2iai2
Vektorové instrukce
podstatná část operací je násobení matic takto implementuje knihovna LAPACK problém algoritmu - „dlouhé" matice A21 a A12 Quintana et. al (2008): algoritmus s bloky 2 p x p
< □ ► 4 s ► 4 i » 4 :
■O0.O
13/18
Paralelní blokové algoritmy
► paralelismus na úrovni BLAS
► existují optimalizované implementace
► implicitní synchronizace na konci každého volání
► nemusí být dosažitelný optimální výkon
► blokové operace jsou efektivní provedené vcelku
► všechna data se dostanou naráz do cache
► na úrovni L1/L2 se vyplatí na jednom jádru CPU
► vzájemně nezávislé blokové operace na více jádrech
Paralelní blokové algoritmy
Quintana et. al (2008) - prototypová implementace SuperMatrix
► konkrétni algoritmus proběhne „abstraktně"
► blokové operace s deklarovanými závislostmi se pouze zařadí do fronty
► následně se vyhodnotí pořadí zpracování
► operace se provedou potenciálně paralelně
čitelná formulace algoritmu bez explicitního paralelismu efektivní provedení
► matice n > 5000
► 16 jader CPU
► LU dekompozice na více než 50% teoretického výkonu
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
15/18
Vektorové instrukce
► historie - např. Cray v předsálí budovy D
► současné vektorové systémy (např. NEC SX)
► řešení velmi specializovaných problémů
► nízký výkon na skalárních operacích
► vektorová rozšíření v architektuře x86
► MMX: pouze celá čísla, kolize s float registry
► SSE: 8 (16) registrů po 128 bitech, postupně přidávané instrukce (rsqrt)
► AVX: 256 bitové registry, 3-operandové instrukce
► snazší využití plného výkonu procesoru
Vektorové instrukce
► triviální příklad
,b[4];
i<4; a[i] += b[i];
32(%rsp), %xmmO (%rsp), %xmmO %xmmO, 32(%rsp)
► vektorizaci kódu provede chytřejší kompilátor (icc)
► musíme hlídat, abychom tomu nebránili
► recyklace proměnných, vedlejší efekty in-line funkcí, ..
► vyplatí se kontrola vygenerovaného assembleru
Využití GPU
téměř vše, co platí o cache, platí také o GPU
► rychlá paměť omezené velikosti
► netriviální režie kopírování z/do hlavní paměti
paralelní kód
► daleko masivnější (desítky až stovky) viz PV197: GPU Programming
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce