vypočtu A. Křenek
PA081: Programování numerických výpočtů
9. Systémy lineárních rovnic a optimalizované implementace lineární algebry
Aleš Křenek jaro 2013
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Motivace
► hledáme řešení systému
a\\X\ + a\zXz
&2\X\ + CI22.X2. &MlXl + ajvf2^2 -
cimnXn = b>m
maticové vyjádření
Ax = b
součást metod řešení nelineárních rovnic, optimalizací atd.
diferenciální rovnice
► simulace dynamických systémů
► metoda konečných prvků
realistické osvětlení scény (radiosita) a mnoho dalších ...
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
■O0.O
Motivace
Simulace pohybu tělesa
► založena na Newtonově zákoně F = ma
► vede na systém diferenciálních rovnic (13 proměnných)
dy(t) dt
( dx/dt \ dq/dt dP/dt y dLjdt i
í   -P \
m \wqq
F T
pro simulaci potřebujeme opakovaně řešit
Yn ~ Yn-1
dt
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
triková aproximace - zanedbání vyšších členů Taylorova rozvoje dy/dt jako funkce y
konečný krok od yn_i k yn znamená řešení systému 13 lineárních rovnic
LU
dekompozice
Vektorové instrukce
Motivace
Modely měkkých tkání
Interakce s měkkými tkáněmi
► tréning chirurgů - vytvořit simulaci chování tkání
► potřebujeme vytvořit matematický model tkáně
► s touto tkání pak můžeme interagovat (např. hapticky)
► je třeba simulovat chování tkáně s dostatečnou přesností
použité matematické modely
► deformovatelná tělesa
► chování je popsáno parciálními diferenciálními rovnicemi, zpravidla neznáme analytické řešení
► řešíme numericky pomocí diskretizace metodou konečných prvků (FEM)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Motivace
Metoda konečných prvků
► simulovaná tkáň je geometricky aproximována meshem
deformace popisuje teorie elasticity
geometricky lineární model - systém lineárních rovnic
geometricky nelineární model - systém nelineárních rovnic
► v každém kroku iterační metody systém lineárních rovnic
systémy rovnic jsou řídké (ovlivňují se jen přímo spojené uzly)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Motivace
Rekonstrukce signálu z ultrazvuku
► pro vývoj filtrů rekonstruujících data z ultrazvuku je výhodné mít počítačový model přístroje a vyšetřovaného objektu
► snáze pak ladíme metody rekonstrukce, nevnášíme nepřesnost danou nedokonalým modelem vyšetřovaného objektu či nedokonalým přístrojem
► pro výpočet šíření vln v objektu použijeme FEM
► na vlnovou délku ultrazvuku potřebujeme několik lineárních elementů
► velikost elementu je pak v desetinách mm
► vyšetřovaný objekt má velikost v jednotkách až desítkách centimetrech
► systémy o jednotkách až desítkách miliónů rovnic
► vysoké nároky na výpočetní kapacitu
► vysoké paměťové nároky
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
□     s   4 :
■OQ.O
Klasifikace problémů
M = N - dobře podmíněné systémy
► naděje na jednoznačné řešení (je-li A regulární) M < N - nedostatečně podmíněné systémy
► žádné řešení
► nekonečně mnoho řešení (N - M-rozměrný prostor) M > N - příliš podmíněné systémy
► přesné řešení zpravidla neexistuje
► hledáme „nejlepší kompromis"
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Klasifikace problémů
M = N - dobře podmíněné systémy
► naděje na jednoznačné řešení (je-li A regulární) M < N - nedostatečně podmíněné systémy
► žádné řešení
► nekonečně mnoho řešení (N - M-rozměrný prostor) M > N - příliš podmíněné systémy
► přesné řešení zpravidla neexistuje
► hledáme „nejlepší kompromis"
husté vs. řídké
► 0(N2) vs. O(N) nenulových prvků
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Klasifikace problémů
M = N - dobře podmíněné systémy
► naděje na jednoznačné řešení (je-li A regulární) M < N - nedostatečně podmíněné systémy
► žádné řešení
► nekonečně mnoho řešení (N - M-rozměrný prostor) M > N - příliš podmíněné systémy
► přesné řešení zpravidla neexistuje
► hledáme „nejlepší kompromis"
husté vs. řídké
► 0(N2) vs. O(N) nenulových prvků
velikost - dopady kumulace chyby
► N ~ 50 - zpravidla stačí f 1 oat
► N -200-500 - double
► větší - vyžadují speciální metody
Přehled metod
prime
► daný algoritmus, vždy stejný počet operací
► nemusí fungovat dobře pro „skoro singulární" nebo příliš velké systémy
► preferované pro „normální" problémy
► Gaussova eliminace, LU dekompozice, ...
iterační
► postupně vylepšované řešení dosažení kritéria konvergence
► různá náročnost pro různé vstupy
► Jacobi, Gauss-Seidel, sdružené směry ...
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Přehled metod
specializované metody pro řídké matice
► pro různé vzory řídkých matic
► výrazně efektivnější, jediná možnost řešení velkých systémů
► přímé i iterační
metody pro singulární systémy
► analyticky i numericky („skoro") singulární
► nalezení celého prostoru řešení
► standardní metody zhavarují
► dekompozice na singulární hodnoty (v příští přednášce) řešení příliš podmíněných systémů
► specializované metody nejmenších čtverců
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Dostupné knihovny
zpravidla sáhneme k hotové knihovně
► nejsme první, kdo tento problém řeší
► implementace optimalizované pro různé případy
► k volbě vhodné metody je třeba rámcově rozumět jejich principům
Numerical Recipes in C
► hlavní zdroj pro tuto přednášku
► jednoduché přehledné implementace
LAPACK http://www.netli b.org/1apack/
► plnohodnotná volně dostupná knihovna
► využívá nízkoúrovňové knihovny BLAS, zpravidla implementované strojově závisle
► obecné a specializované funkce pro různé typické případy
NAG http://www.nag.co.uk/
► rozsáhlá komerční numerická knihovna a další...
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Gaussova eliminace
Základní postup
standardní „školní" metoda
► řešení systému se nezmění, nahradíme-li libovolný řádek A a odpovídající prvek b lineární kombinací tohoto řádku
s libovolným dalším
vynulujeme 0,21, ■ ■ ■, &m\ odečtením ^ násobku prvního řádku
obdobně pokračujeme druhým sloupcem od a^2
výsledkem je matice s vynulovanými prvky pod diagonálou
řešení systému získáme zpětnou substitucí
► poslední rovnice umm^m =     je triviální
► do předposlední dosadíme získanou hodnotu xm atd.
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
□ 4 :
Gaussova eliminace
Numerické problémy
► diagonální prvek atk je 0 nebo příliš malý
► v k-té iteraci se jím dělí
► dojde k přetečení
► při odečítání dochází k přílišné ztrátě platných cifer
6 1 1 1
1
0 1
zaokrouhlení může vést až k azz = -7, to odpovídá původní matici
" e 1 1 0
metoda v této podobě je numericky nestabilní
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Gaussova eliminace
Pivotíng
► další tvrzení o systému
► řešení systému se nezmění výměnou libovolné dvojice řádků
► řešení systému se nezmění výměnou libovolné dvojice sloupců, zaměníme-li zároveň příslušné komponenty vektoru x
► pivot je prvek, který po aplikaci těchto kroků použijeme k dělení při eliminaci fe-tého sloupce
► částečný pivotíng (parciální)
► v iteraci k vybíráme z ojt pro j > k, tj. jen z nulovaného sloupce
► plný pivoting
► vybíráme z ay pro i, j > k, tj. z celé podmatice doprava a dolů
► vyžaduje udržovat vznikající permutaci proměnných
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Gaussova eliminace
Pivotíng
► za pivota volíme maximum z kandidátů
► pro všechny faktory v elirninačních krocích platí | ^ | < 1 parciální i plný pivotíng jsou pak numericky stabilní přínos plného pivotingu je jen okrajový
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Gaussova eliminace
Pivoting
► za pivota volíme maximum z kandidátů
► pro všechny faktory v eliminačních krocích platí | ^ | < 1
► parciální i plný pivoting jsou pak numericky stabilní
► přínos plného pivotingu je jen okrajový
► volba pivota závisí na řádu koeficientů původního systému
► vynásobení jedné rovnice faktorem 1010 ...
► za pivota lze volit koeficient, který by byl největší, kdyby byl celý původní systém normalizovaný, tj. největší koeficient všech rovnic roven 1
► vyžaduje dodatečnou údržbu faktorů škálování, může přinést větší robustnost
Gaussova eliminace
Gauss-Jordánova eliminace
Gaussovu metodu lze použít pro více pravých stran současně
speciálně pro N bázových vektorů dostaneme výpočet matice A-1
rozšíření - Gauss-Jordanova eliminace
► v každé iteraci eliminujeme prvky nad i pod diagonálou A
► výsledkem je jednotková matice
► řešení systému je přímo modifikovaný vektor b
► resp. dostáváme A-1
srovnatelné se základní implementací
► reálně provádíme tytéž operace
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Gaussova eliminace
Výhody a nevýhody
► jednoduchá, dobře pochopitelná a stabilní metoda
►•je třeba znát pravé strany dopředu
► jsou součástí celého výpočtu
► výpočet A-1 je zatížen poměrně velkou numerickou chybou
► přímý výpočet řešení systému Ax = b je přesnější než výpočet A-1 a následné x = A_1b
Rozklady matic
rozklad matice na vhodné činitele
► řešení systému rovnic je pak triviální LU
► čtvercové matice
► de facto jiné vyjádření Gaussovy eliminace
► explicitní vyjádření při více pravých stranách
Choleského
čtvercové, symetrické, pozitivně definitní matice
► stabilnější a rychlejší algoritmus
QR
► obecné mxn
► numerický stabilní, náročnější výpočet singulární hodnoty
► náročnější výpočet, extrémně stabilní
► dává přímo „řešení" pro M $ N
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Iterační zpřesnění
► i přímé metody řešení lineárních rovnic ztrácí přesnost
► v závislosti na velikosti systému a poměru koeficientů
► kumulace chyb pocházejících ze sčítání/odčítání
► v optimistickém případě 2-3 platná místa
► u „skoro singulárních" matic podstatně horší
► lze vylepšit iteračním procesem
► x je přesné řešení Ax = b, známe ale jen x + <5x; položíme
A(x + 5x) = b + 5b
odečtením dostaneme
A5x = 5b = A(x + 5x) - b
pravou stranu umíme spočítat, řešíme nový systém rovnic pro 5x
► pro výpočet pravé strany je nutná vyšší přesnost odčítání
► s výhodou využijeme recyklaci LU dekompozice
výsledným 5x korigujeme původní x
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Iterační zpřesnění
► schematické znázornění
-»-f A
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
postup lze opakovat, zpravidla stačí 1-2 krát
LU
dekompozice
Vektorové instrukce
Řídké matice
typicky rozsáhlé problémy (miliony rovnic)
► ale počet nenulových koeficientů je malý, zpravidla O(N) v extrémním případě neřešitelné standardními metodami
► příliš velká paměťová a časová náročnost
► de-facto nepřekonatelné problémy s přesností výpočtu
specializované metody
► využívají nulových koeficientů
► vyžadují jen O (N) operací i paměti
► závislé na konkrétním vzoru nenulových prvků
► např. LU dekompozice tridiagonální matice - jen N prvkový vektor navíc a 2 cykly o N - 1 iteracích
v některých případech aproximační
► nenulové prvky se během řešení „množí" nad míru
► je třeba některé zanedbávat
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Řídké matice
Některé speciální vzory
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Řídké matice
Uložení v paměti
problematická není jen výpočetní náročnost, ale zejména nároky na paměť
► N = 106 by vyžadovalo ve floatech 4 TB existují různá schémata, např.
► pole hodnot f 1 oat val [] a indexů i nt idx[]
► prvních N prvků val [] jsou všechny diagonální prvky, zpravidla bývají nenulové
► prvních N prvků i dx [] jsou indexy ve val [], kde jsou uloženy první nenulové nediagonální prvky jednotlivých řádků matice
► i dx [N + 1] ukazuje za poslední prvek val []
► další prvky val [] jsou nenulové nediagonální prvky matice uspořádané po řádcích a sloupcích
► další prvky i dx [] jsou indexy sloupců odpovídajících prvků val []
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
□     s   4 :
■OQ.O
Řídké matice
Uložení v paměti
původní matice
je reprezentována
3	0	1	0	0
0	4	0	0	0
0	7	5	9	0
0	0	0	0	2
0	0	0	6	5
	1	2	3	4	5	6	7	8	9	10	11
idx[]	7	8	8	10	11	12	3	2	4	5	4
val[]	3	4	5	0	5	n/a	1	7	9	2	6
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Řídké matice
Sherman-Morrisonův postup
► řídké matice, které se „trochu liší" od známého vzoru
► navíc řádek, sloupec apod.
► obecně A' = A + (u ® v)
► lze ukázat
(A'
(A + (u ® v)
(A-1!!) 0 (vA'1) 1 + vA_1u
některou z hotových metod vypočteme A-1 aplikací vzorce získáme (A')-1
► je-li A-1 řídká v řádu O(N), náročnost celé metody je O(N) postup lze opakovat pro různá u, v existují další zobecnění (Woodbury, viz literatura)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Iterační metody
► obecně méně přesné než přímé metody
► řešení řídkých systémů
► neexistuje přímá metoda pro konkrétní vzor
► výpočet iteračního kroju je zpravidla O(N)
► nevyžaduje další paměť
► téměř singulární systémy
► přímé metody ztrácí přesnost kumulací chyby
Iterační metody
Prostá iterace
rovnici Ax = b převedeme na tvar x = A'x + b' opakovaně počítáme iterační krok kritérium konvergence
► zobrazení x <- A'x + b' musí být kontrakce
► stejné jako u nelineárních rovnic
naplnění předpokladů věty o pevném bodě
► lmifc^K. l(A')fc| = 0 pro nějakou maticovou normu
► Frobeniova norma _
|A|
spektrální norma
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
|A| =p(A7A)
kde pO je maximální absolutní hotnota vlastních hodnot matice
maximální součet sloupce nebo řádku
LU
dekompozice
Vektorové instrukce
Iterační metody
Prostá iterace - příklady konkrétních metod
Jacobi: z rozkladu A = D - L - U dostaneme x = D X(L + U)x + D xb, tj.
Gauss-Seidel: x = (D - L)_1Ux + (D - L) ^Ub
,    / i—l n
au \      j=i j=í+i
lze recyklovat uložení x
konvergence není zaručena automaticky
► pro konkrétní A některé metody konvergují, jiné ne
► specifická kritéria viz literatura
Iterační metody
Metoda konjugovaných gradientů
► lze využít i metody řešení nelineárních rovnic
► smysl to má jen ve velmi specifických případech
► výjimkou je metoda konjugovaných gradientů
► minimalizujeme funkci
/(x)
1
xAx - bx
v minimu je její gradient nulový, tj.
0 = V/ = Ax - b
tedy minimum / přesně odpovídá řešení Ax = b metoda tedy najde minimum po N iteracích
► / je kvadratická forma, viz minulá přednáška
takto jednoduše funguje jen pro pozitivně definitní A, rozšířit
při výpočtu je třeba jen násobit Ax
► to lze u řídkých matic velmi efektivně
lze
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Lineární algebra rychle
proč má smysl optimalizovat právě algoritmy lineární algebry?
frekventované problémy
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
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
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]
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
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. ...
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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 i = y i + a
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 vektory spojitě, ne jako sloupce matice (v C)
► jinak se s tím nedá nic moc dělat
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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í
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Asymptotická složitost
► časová složitost násobení matic je 0{n3)
4 a ►   4 (5 ►   <       ►   4 S ►
PA081: Programování numerických
výpočtů
A. Křenek
Motivace
Přehled metod
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
B LAS
Asymptotická složitost
LU
dekompozice
Vektorové instrukce
•o<\o
39/46
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 + 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)
0(n
2.81 >
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
□ 4 :
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 + 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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Domácí úkol
implementujte co nejefektivněji Strassenův algoritmus
od určité velikosti matice přepněte na BLAS
vyhodnoťte chování (kdy začne být výhodnější, co to stojí)
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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 :
LU
dekompozice
Vektorové instrukce
■O0.O
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 problém algoritmu - „dlouhé" matice A21 a A12 Quintana et. al (2008): algoritmus s bloky 2 p x p
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
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
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce
Shrnutí
řešení lineárních rovnic je součást řady problémů
v operacích LA se tráví velká část výpočetního času
má smysl hledat optimalizované a numericky stabilní algoritmy
neexistuje univerzální řešení
k dispozici řada optimalizovaných implementací
Gaussova eliminace
Iterační zpřesnění
Řídké matice
Iterační metody
LA rychle
Vyrovnávací paměti
Blokové algoritmy
LU
dekompozice
Vektorové instrukce