PA081: Programování numerických
výpočtů A. Křenek
PA081: Programování numerických výpočtů
8. Lineární algebra rychle
Aleš Křenek
jaro 2011
1/18
Motivace
Vyrovnávací paměti
Blokové algoritmy
BLAS
Asymptotická složitost
LU
dekompozice
Vektorové instrukce
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ů
PA081: Programování numerických
výpočtů
A. Křenek
Vyrovnávací paměti
Blokové algoritmy
Asymptotická složitost
LU
dekompozice
Vektorové instrukce
2/18
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
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
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
double pol e[100] [512][8] ; f or (i=0; i<100; i++) a += pol e [i ] [0] [0] ;
► řešením je deklarace pol e [100] [513] [8]
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
PA081: Programování numerických
výpočtů
A. Křenek
Motivace
Vyrovnávací paměti
Blokové algoritmy
BLAS
Asymptotická složitost
LU
dekompozice
Vektorové instrukce
5/18
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 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
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 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
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 f or í = 1.....n
načti řádek Aj* do cache
f o r j = 1.....n
yi = yi + AijXj
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/m ~ 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í matíc C = C + AB
f or i = 1.....n
načti řádek Aj* do cache f o r / = 1.....n
načti Qj do cache
načti sloupec B* j do cache for k = 1.....n
= Qj + AíkBkj
zapiš Qj do hlavní paměti
► přístupy do paměti m = n2 + n3 + 2n2
► aritmetické operace / = 2n3
► efektivita opět « 2
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 x b, b = n/N
do cache
f or i = 1.....JV
f or j = 1.....N
načti blok C for k = l,..'.,N
načti blok Ajfc do cache načti blok By do cache Cjj = Cjj + AjfcBfcj (maticové násobení) zapiš blok Cjj do cache
► přístupy do paměti
► čtení bloků A: N3 (n/N)2 = N * n2 - dtto B: N * n2
► čtení a zápis C: 2n2
2n3
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
efektivita f /m
(2JV+2)n2
b
9/18
Blokové algoritmy
Maticové násobení
► velikost reálne Ll cache je 32 kB,
► do Ll se vejdou 3 matice velikosti 36 x 36 (v doubl e)
► 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á
► žá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í
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
10/18
Knihovna BIAS
► Basic Linear Algebra Subprograms, http://www.netli b.org/bl as
► základní operace, např.
► DOT: xy
► AXPY: y + Ax
► NRM2: |x|2
► TRMV: Ax
► TRSV: A_1x
► CEMM: PC + aAB
► verze pro typy f 1 oat,doubl e,compl ex
► specializované varianty pro symetrické, trojúhelníkové, 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)
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
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)(Bn	+ B22)	
M2	= (A2i + A22)Bn		
			C11
M3	= An(Bi2 + B22)		
			C12
M4	= A22(B2i -Bn)		
			C21
M5	= (Au + Ai2)B22		
			C22
M6	= (A2i -AiiKBn	+ Bi2)	
M7	= (A12 - A22)(B2i	+ B22)	
Mi + M4 - M5 + M7 M3 + M5 M2 + M4
Mi - M2 + M3 + M6
rekurzivní aplikace, počet operací T(n)
T(n) = 7T(n/2) + 18(n/2)2 = 0(nlo&7) = 0(n281)
A. Křenek
Motivace
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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)(Bn	+ B22)	
M2	= (A2i + A22)Bn		
			C11
M3	= An(Bi2 + B22)		
			C12
M4	= A22(B2i -Bn)		
			C21
M5	= (Au + Ai2)B22		
			C22
M6	= (A2i -AiiKBn	+ Bi2)	
M7	= (A12 - A22)(B2i	+ B22)	
Mi + M4 - M5 + M7 M3 + M5 M2 + M4
Mi - M2 + M3 + M6
rekurzivní aplikace, počet operací T(n)
T(n) = 7T(n/2) + 18(n/2)2 = 0{nXo^7) = 0(n2-81) ► Coppersmith-Winograd (1987): 0(n2376)
A. Křenek
Motivace
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
LU dekompozice
► rekurzivní formulace algoritmu
/A
oo
aío
aoi   A02 \
«11 a2i
\ A20
skalární a bloková verze
An
.A2i A12 A22
a12
A22 /
= a2i/«n a{2 zůstává A22 = A22 - a2iai2
a2i
LU
A11 A21 Lx/Ais
A22 - L21A12
► podstatná část operací je násobení matic
► takto implementuje knihovna LAPACK
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
13/18
LU dekompozice
► rekurzivní formulace algoritmu
/A
oo
aío
aoi   A02 \
«11 a2i
\ A20
► skalární a bloková verze
An
.A2i A12 A22
a12
A22 /
= a2i/«n a{2 zůstává A22 = A22 - a2iai2
a2i
LU
A11 A21 Lx/Ais
A22 - L21A12
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 2p x p
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
PA081: Programování numerických
výpočtů
A. Křenek
Motivace
Vyrovnávací paměti
Blokové algoritmy
BLAS
Asymptotická složitost
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
► 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
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
17/18
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