PA081: Programování numerických
výpočtů
A. Křenek
PA081: Programování numerických výpočtů
Aleš Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
jaro 2011
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
1/38
O čem to bude
► zajímá nás chování skutečného světa
► problémy přírodovědné, technické, humanitní...
► popisujeme matematickými prostředky
► zejména pomocí reálných čísel
► umělý aparát, leckdy zcela neodpovídá skutečnosti
► za staletí celkem zvládnutý, vyučovaný od základní školy, obecně přijímaný
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
2/38
O čem to bude
► zákonitosti zkoumaných systémů typicky vyjádřeny rovnicemi
► zajímavé systémy ->• složité rovnice
► reprezentativní příklad - Schrodingerova rovnice
íh^-Y(x,t) =ÍT¥(x,t) ot
► analyticky řešitelná jen pro triviální případy
► izolovaná částice, částice v potenciálové jámě.....atom
vodíku
► i tak je řešení dost komplikované
► numerické řešení je jediný realisticky možný přístup
O čem to bude
► numerická matematika
► řešení matematicky formulovaného problému aritmetickými prostředky
► zpravidla algoritmický postup - numerická metoda
► disciplina podstatně starší než počítače
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
4/38
O čem to bude
► numerická matematika
► řešení matematicky formulovaného problému aritmetickými prostředky
► zpravidla algoritmický postup - numerická metoda
► disciplina podstatně starší než počítače
metody jsou známé, naprogramujeme je, a je hotovo
► ne tak docela
(De)motivační příklad
PA081: Programování numerických
výpočtů
A. Křenek
► pro dané n vypočítejte integrál
O čem to bude
(De)motivační příklad
ri
Obsah předmětu
xnex-l
dx
Předpoklady návaznosti
Jo
Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
5/38
(De)motivační příklad
► pro dané n vypočítejte integrál
-i
xnex ldx
o
► očekávané vlastnosti
► E0 = [e^}\ = l-\
► pro 0 < x < 1 platí xn > 0 a 0 < ex~l < 1, tedy En > 0
podobně
r1 1
£n <    xndx =-- a tedy   lim En = 0
Jo n + 1
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
5/38
(De)motivační příklad
► integrací per-partes u = xn, dv = ex ldx
-i
En = [xnex-l]\
nxn xex 1dx = 1 - nE^-
o
n-l
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
6/38
(De)motivační příklad
► integrací per-partes u = xn, dv = ex ldx
-i
En = [xnex-l]\
nxn xex 1dx = 1 - nE^-
o
n-l
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
6/38
(De)motivační příklad
1: 0.367879
2: 0.264241
3: 0.207277
4: 0.170893
5: 0.145534
6: 0.126796
7: 0.112430
8: 0.100563
9: 0.094933
10: 0.050674
11: 0.442581
12: -4.310974
13: 57.042664
14: -797.597290
15: 11964.958984
16: -191438.343750
17: 3254452.750000
18: -58580148.000000
19: 1113022848.000000
(De)motivační příklad
► co je špatné?
► pracujeme s konečnou reprezentací reálného čísla
► společný problém ručního i strojového zpracování
► v počítači daleko plíživější podoba (nevidíme mezivýsledky)
výpočtu A. Křenek
(De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
8/38
(De)motivační příklad
► co je špatně?
► pracujeme s konečnou reprezentací reálného čísla
► společný problém ručního i strojového zpracování
► v počítači daleko plíživější podoba (nevidíme mezivýsledky)
► v proměnné E [n] není uloženo přesně En
už na začátku zatíženo chybou, E[0] = Eq + e
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
8/38
(De)motivační příklad
► co je špatně?
► pracujeme s konečnou reprezentací reálného čísla
► společný problém ručního i strojového zpracování
► v počítači daleko plíživější podoba (nevidíme mezivýsledky)
► v proměnné E [n] není uloženo přesně En
už na začátku zatíženo chybou, E[0] = Eo + e
► i při zcela přesném výpočtu:
E [1] = 1 - E [0] = 1 - Eo - e = Ei - e
E[2] = 1 - 2E[1] = 1 - 2Ei + 2e = E2+2e
E[3] = 1 - 3E[2] = 1 - 3E2 - 3(2e) = £3 - 3(2e)
E[n] =
= En + (-l)nn\e
O čem to tedy bude
► představení numerických metod pro řešení vybraných problémů
► pragmaticky, bez rozsáhlé teoretické analýzy, okrajových podmínek atd.
► formulace přiměřeně přesného a efektivního algoritmu
► matematicky korektní metoda nestačí
► pro různá použití jsou vhodné různé alternativy
► použití na smysluplném příkladu
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
9/38
Obsah předmětu
► obecné zásady psaní numerického kódu
► nelineárních rovnice
f (x) = y
► optimalizace (hledání lokálního minima)
► lineární úlohy
Ax = B
► automatické derivování
► diferenciální rovnice (zlehka)
Předpoklady a návaznosti
► předpokládané znalosti
► návrh algoritmů, elementární schopnost programování
► porozumět kódu v C
► programovat ve svém oblíbeném jazyce (C++, Java, Fortran, Python, ...)
► základy lineární algebry a matematické analýzy
► základní orientace v architektuře a assembleru x86
► http://www.i ntel.com/products/processor/manuals/
► pro referenci, netřeba číst všechno ;-)
Předpoklady a návaznosti
► předpokládané znalosti
► návrh algoritmů, elementární schopnost programování
► porozumět kódu v C
► programovat ve svém oblíbeném jazyce (C++, Java, Fortran, Python, ...)
► základy lineární algebry a matematické analýzy
► základní orientace v architektuře a assembleru x86
► http://www.i ntel.com/products/processor/manuals/
► pro referenci, netřeba číst všechno ;-)
► návaznosti
► M4180: Numerické metody I
► PV027: Optimalizace
► IA039: Architektura superpočítačů a intenzivní výpočty
► PV192: Paralelní technické systémy
► PV197: GPU Programming
► diplomové práce
Zkouška
► znalosti v rozsahu přednášek
► písemná část
► teorie i praktický příklad
► záměrně tvrdší hodnocení (E > 50%)
► ústní část
► šance na vylepšení
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
12/38
Domácí úkoly a konzultace
► domácí úkoly
► dobrovolné, slouží k lepšímu pochopení problema
► zadávány průběžně několikrát za semestr
► řešení do dalšího týdne
► drobné body k dobru ke zkoušce
► konzultace
► v mezích možností kdykoli po domluvě emailem
► ÚVT, Botanická 68a, kancelář C122
Čísla v plovoucí řádové čárce
► standard IEEE 754
vychází návrhu reprezentace čísel v koprocesoru Intel 8087
► základní formát ±l.mmmmmmmmm... x 2±ee-
► znaménko mantisy (1 bit)
► mantisa l.mmmmmmmm
► binárně, absolutní hodnota
► číslo v intervalu [1, 2) v dané přesnosti
► nekonečná množina pokryta konečným počtem hodnot
► základní zdroj nepřesnosti
► exponent
► binární číslo v rozsahu 1 až např. 254
► znamená hodnoty exponentu -126 až +127
► speciální význam hodnot 00... 00 a 11... 11 - kódování ±0, ±oo, ±NaN
Čísla v plovoucí řádové čárce
► nejběžnější typy - velikosti v bitech a rozsah exponentu
typ	exp.	man tisa	celkem	rozsah exp.
Single (float)	8	23	32	127
Double	11	52	64	1023
Quad (long double)	15	112	128	16383
► tj. např. největší číslo typu float je 2127 = 10
► nejmenší nenulové 2~126 = 10~37
► relativní chyba je omezena 2~25 = 10~8
► tedy počítáme na cca. 8 platných cifer
Problémy konečné reprezentace
Ztráta presnosti sčítání
► pro zjednodušení v desítkové soustavě a na 3 platné cifry man tisy
► sečtěme 1.23e3 a l.OOeO
► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3
► v rámci dané přesnosti korektní výsledek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
16/38
Problémy konečné reprezentace
Ztráta přesnosti sčítání
► pro zjednodušení v desítkové soustavě a na 3 platné cifry man tisy
► sečtěme 1.23e3 a l.OOeO
► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3
► v rámci dané přesnosti korektní výsledek
► k 1.23e3 desetkrát přičtěme l.OOeO
► opakované aplikování předchozího postupu dává opět 1.23e3
► nemusí to být to, co jsme chtěli
► 1.23e3 + (l.OOeO + ■ ■ ■ + l.OOeO) = 1.24e3
Problémy konečné reprezentace
Ztráta přesnosti sčítání
► pro zjednodušení v desítkové soustavě a na 3 platné cifry man tisy
► sečtěme 1.23e3 a l.OOeO
► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3
► v rámci dané přesnosti korektní výsledek
► k 1.23e3 desetkrát přičtěme l.OOeO
► opakované aplikování předchozího postupu dává opět 1.23e3
► nemusí to být to, co jsme chtěli
► 1.23e3 + (l.OOeO + ■ ■ ■ + l.OOeO) = 1.24e3
► operace v plovoucí řádové čárce nejsou asociativní
► ani tam, kde jsme na to z algebry zvyklí
Problémy konečné reprezentace
Ztráta presnosti sčítání
O čem to bude (De)motivační príklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Dobrá rada #1
Čekáme-li kumulativní efekt, sčítání a odčítání je třeba provést nejprve na číslech se srovnatelným exponentem.
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
PA081: Programování numerických výpočtů
A. Křenek
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
17/38
Problémy konečné reprezentace
Katastrofální zrušení
► odečtení dvou vzájemně blízkých čísel vede k výrazné ztrátě přesnosti
► počítáme např. 1.234e3 - 1.218e3
► na 3 cifry reprezentováno jako 1.23e3 a 1.22e3
► relativní chyba 0.32% a 0.16%
► odečtením dostaneme l.OOel
► správný výsledek byl 1.600el, relativní chyba 38%
► navíc působí dojmem falešné přesnosti
Problémy konečné reprezentace
Katastrofální zrušení
► obráceně
► z předchozího výpočtu máme výsledky na 3 místa: 1.23e3 a 1.22e3
► takto je vypisuje pri ntf ("%8. 2e" , . . . )
► v proměnných ale je 1.234e3 a 1.218e3
► očekávaný výsledek odečtení: lei
► dostaneme ale 1.6el
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
19/38
Problémy konečné reprezentace
Katastrofální zrušení
► obráceně
► z předchozího výpočtu máme výsledky na 3 místa: 1.23e3 a 1.22e3
► takto je vypisuje pri ntf ("968. 2e" , . . . )
► v proměnných ale je 1.234e3 a 1.218e3
► očekávaný výsledek odečtení: lei
► dostaneme ale 1.6el
Dobrá rada #2
Vyhněme se odčítání dvou blízkých čísel. Je-li to i tak nezbytné, počítejme s výsledkem zatíženým potenciálně velkou chybou.
Problémy konečné reprezentace
Násobení a dělení
► samo o sobě nezanáší významnou chybu
► zachovává existující chybu
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
20/38
Problémy konečné reprezentace
Násobení a dělení
► samo o sobě nezanáší významnou chybu
► zachovává existující chybu
► příklad: (a + b)2 pro a = 1.23, b = 0.0155, na 3 cifry přesný výsledek je 1.55127025
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
20/38
Problémy konečné reprezentace
Násobení a dělení
► samo o sobě nezanáší významnou chybu
► zachovává existující chybu
► příklad: (a + b)2 pro a = 1.23, b = 0.0155, na 3 cifry přesný výsledek je 1.55127025
► přímý výpočet na 3 platné cifry:
a + b = 1.23 + 0.02 = 1.25
► tedy (a + b)2 = 1.56, chyba 0.56%
► není to tak zlé, ale může být lepší
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
20/38
Problémy konečné reprezentace
Násobení a dělení
► po transformaci (a + b)2 = a2 + 2ab + b2:
1.51 + 2 X 0.0191 + 0.000240 =
1.51 + 0.0382 + 0.000240 = 1.55
► chyba 0.082%, tedy téměř o řád menší
► za cenu 6 aritmetických operací místo 2
► bude 3x pomalejší
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
21/38
Problémy konečné reprezentace
Násobení a dělení
► po transformaci (a + b)2 = a2 + 2ab + b2:
1.51 + 2 x 0.0191 + 0.000240 =
1.51 + 0.0382 + 0.000240 = 1.55
► chyba 0.082%, tedy téměř o řád menší
► za cenu 6 aritmetických operací místo 2
► bude 3x pomalejší.. .uvidíme
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
21/38
Měření rychlosti výpočtu
Naivní přístup
► POSIX volání gettimeofdayO
gettimeofday(&start,NULL); c = a + b; c *= c;
getti meofday(&stop,NULL);
► současné CPU 2-3 GHz
► 1 aritmetická operace ~ l-10ns nemáme tak přesné hodiny
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
22/38
Měření rychlosti výpočtu
Zopakujeme v cyklu
► opakování 109x navýší čas na měřitelnou
gettimeofday(&start,NULL); for (i=0; i<1000000000;
c = a + b;
c *= c;
}
getti meofday(&stop,NULL);
Měření rychlosti výpočtu
Zopakujeme v cyklu
► opakování 109x navýší čas na měřitelnou
gettimeofday(&start,NULL); for (i=0; i<1000000000; i++)
c = a + b;
c *= c;
}
getti meofday(&stop,NULL);
počítá rychlostí 18 Gflop/s
► měření na Intel Xeon E5520 2.27 GHz
► trochu podezřelé
Měření rychlosti výpočtu
Zopakujeme v cyklu
► optimalizující kompilátor gcc -03 -f unrol 1-1 oops
movss	44(%rsp), %xmmO
xorl	%eax, %eax
addss	40(%rsp), %xmmO
mul ss	%xmmO, %xmmO
addl	$8, %eax
cmpl	$1000000000, %eax
j ne	.L2
movss	%xmm0, 36(%rsp)
v těle cyklu se nic nepočítá optimalizaci obecně chceme
► použití registrů pro proměnné, max. využití FPU, eliminace zbytečného opakování je příliš
►  4 J ►
Měření rychlosti výpočtu
Zopakujeme v cyklu
► brzdící kód
► uměle vložený do těla cyklu
► za běhu programu se nevyvolá - nebrzdí doopravdy
► zabrání příliš agresivní optimalizaci cyklu
► kompilátor musí předpokládat:
► proměnná ni kdy nemusí být 0
► funkce brzdaQ má vliv na odkazované proměnné
Měření rychlosti výpočtu
Jak to dopadlo
► Intel XeonE5520 2.27GHz
vzorec	čas	Mflop/s	cyklů/s
(a + b)2	1.33	1509	754
a2 + 2ab + b2	1.99	3018	503
► rozvinutý výpočet jen 1.5x
► vnitřní paralelismus procesoru
► více FPU jednotek
► pipelining
► spekulativní výpočet větví programu
Problémy konečné reprezentace
Násobení a dělení
Dobrá rada #3
Snažme se vzorce transformovat tak, aby se nejdříve násobilo/dělilo, pak teprve sčítalo/odčítalo. Je šance na přesnější výsledek, dopad na výkon nemusí být významný.
Domácí úkol
Funkce brzda() v předchozím příkladu pracuje i s proměnnou ni kdy.
Změňte jen na brzda(&a,&b,&c), vyzkoušejte, co se stane, a vysvětlete.
Bude třeba podívat se do vygenerovaného assembleru.
Původní zdrojový text ve studijních materiálech na ISu, adresář /um/pri klady.
Odevzdání do 2.2.2011, 2 body (96)
Problémy konečné reprezentace
Přetečení a podtečení
► násobení dvou velkých čísel může vést k nereprezentovatelnému exponentu
1.2e30 x 1.2e30 = 1.44e61
► typ fl oat umí exponent [-37, 38]
► výsledek operace už je oo
Problémy konečné reprezentace
Přetečení a podtečení
► násobení dvou velkých čísel může vést k nereprezentovatelnému exponentu
1.2e30 x 1.2e30 = 1.44e61
► typ f 1 oat umí exponent [-37, 38]
► výsledek operace už je oo
► podobně násobení dvou malých čísel vede k podt
► podle nastavení aritmetiky je výsledek ±0 nebo denormalizované číslo
O.OOOOOOmmmx 1(T37
► další výpočty nepřesné a navíc citelně pomalejší
Problémy konečné reprezentace
Přetečení a podtečení
Dobrá rada #4
Při násobení a dělení více čísel uspořádejme posloupnost operací, abychom nedělili velmi malé číslo velkým, velké malým, a nenásobili vzájemně velká nebo malá čísla.
Problémy konečné reprezentace
Porovnávání
► nevinný výpočet
float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2;
printf("%f %f %s\n",c,d,c == d ?
"jsou stejná" : "jsou ruzna");
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
31/38
Problémy konečné reprezentace
Porovnávání
► nevinný výpočet
float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2;
printf("%f %f %s\n",c,d,c == d ?
"jsou stejná" : "jsou ruzna");
má překvapivý výstup
3.010000 3.010000 jsou ruzna
Problémy konečné reprezentace
Porovnávání
nevinný výpočet
float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2;
printf("%f %f %s\n",c,d,c == d ?
"jsou stejná" : "jsou ruzna");
► má překvapivý výstup
3.010000 3.010000 jsou ruzna
► přesnější formát výpisu "%15.12f":
3.009999990463 3.010000228882 jsou ruzna
Problémy konečné reprezentace
Porovnávání
► operátor == na reálných číslech nemá valný smysl
► téměř vždy je zasažen chybou předchozího výpočtu
► místo něj test na dostatečnou blízkost
fabs(c-d) < EPSILON
► stanovení toleranční konstanty zpravidla empiricky
► různá pro hodnoty le-15 a le30
► silně záleží na dané úloze
► ovlivněna i každým konkrétním výpočtem
Problémy konečné reprezentace
Porovnávání
► porovnávání < a > trpí stejným problémem
► záludnější podoba a komplikovanější řešení
► chci porovnat přesné hodnoty x a y
► znám jen přibližná (spočtená) x = x ± ex, y = y ± ey:
i-•-1 x
X
y
y
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
33/38
Problémy konečné reprezentace
Porovnávání
► porovnávání < a > trpí stejným problémem
► záludnější podoba a komplikovanější řešení
► chci porovnat přesné hodnoty x a y
► znám jen přibližná (spočtená) x = x ± ex, y = y ± ey\
"-•-" x
X
y
y
> opatrný („miserly") přístup
► intervaly se i částečně překrývají => tvrdíme x = y
► jinak porovnáme x $ y
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Katastrofální zrušení
Měření rychlosti výpočtu
Přetečení a podtečení Porovnávání
Numerická stabilita
33/38
Problémy konečné reprezentace
Porovnávání
► porovnávání < a > trpí stejným problémem
► záludnější podoba a komplikovanější řešení
► chci porovnat přesné hodnoty x a y
► znám jen přibližná (spočtená) x = x ± ex, y = y ± ey:
i-•-1 x
x
I-•-1 y
y
► opatrný („miserly") přístup
► intervaly se i částečně překrývají => tvrdíme x = y
► jinak porovnáme x $ y
► hladový („eager") přístup
► chceme vědět, že mohlo nastat x > y
► porovnáváme x + ex > y - ey.
Problémy konečné reprezentace
Porovnávání
Dobrá rada #5
Na rovnost porovnávejme vždy s tolerancí.
U porovnaní na nerovnost si uvědomme, co chceme vědět.
Numerická stabilita
(Pseudo)definice
► „metoda počítá sice špatně, ale jenom trochu špatně"
► žádoucí a přitom realistická vlastnost všech numerických algoritmů
Numerická stabilita
(Pseudo)definice
► „metoda počítá sice špatně, ale jenom trochu špatně"
► žádoucí a přitom realistická vlastnost všech numerických algoritmů
► pro vstup x hledáme řešení y = f (x)
► ve skutečnosti
► na aproximaci vstupu x = x + ex
► nepřesnou numerickou metodou /
► dostaneme výsledek y = y + ey
x
Numerická stabilita
(Pseudo)deflnice
► „metoda počítá sice špatně, ale jenom trochu špatně"
► žádoucí a přitom realistická vlastnost všech numerických algoritmů
► pro vstup x hledáme řešení y = f (x)
► ve skutečnosti
► na aproximaci vstupu x = x + ex
► nepřesnou numerickou metodou /
► dostaneme výsledek y = y + ey
► algoritmus je stabilní, ^ existuj e-li pro každé x malé ex \ / tak, že ey je malé x
Numerická stabilita
(Pseudo)definice
► „metoda počítá sice špatně, ale jenom trochu špatně"
► žádoucí a přitom realistická vlastnost všech numerických algoritmů
► pro vstup x hledáme řešení y = f (x)
► ve skutečnosti
► na aproximaci vstupu x = x + ex
► nepřesnou numerickou metodou /
► dostaneme výsledek y = y + ey
► algoritmus je stabilní, existuj e-li pro každé x malé ex tak, že ey je malé
silnější definice zpětné stability - ey = 0
Numerická stabilita
(Pseudo)definice
► „metoda počítá sice špatně, ale jenom trochu špatně"
► žádoucí a přitom realistická vlastnost všech numerických algoritmů
► pro vstup x hledáme řešení y = f (x)
► ve skutečnosti
► na aproximaci vstupu x = x + ex
► nepřesnou numerickou metodou /
► dostaneme výsledek y = y + ey
► algoritmus je stabilní,
existuj e-li pro každé x malé ex tak, že ey je malé
silnější definice zpětné stability
~-y
0
► pro speciální druhy problémů jiné definice stability ► např. numerické integrování - „systém nevybuchne"
Numerická stabilita
Příklady
► addss je numericky stabilní implementace sčítání
► relativní chyba 5 sčítání v typu f 1 oat je omezena
\5\ < 2~25
y = adds(%i, %2) = {x\ +X2)(1 + <5) = x\(l + 5) +X2Í1 + tj. našli jsme
e*i,2 = xi,2<5  a  ey = (xi + x2)5 vyhovující definici
Numerická stabilita
Príklady
► sčítání řady čísel 12 30 + 1 + 1 + 1 + ...je nestabilní
► stabilní alternativa - setřídit a začít od nejmenšího
PA081: Programování numerických
výpočtů
A. Křenek
O čem to bude (De)motivační příklad
Obsah předmětu Předpoklady a návaznosti Zkouška
Čísla v plovoucí řádové čárce
Reprezentace
IEEE 754
Nepřesné sčítání
Katastrofální
zrušení
Násobení a dělení Měření rychlosti výpočtu Domácí úkol Přetečení a podtečení Porovnávání
Numerická stabilita
Shrnutí
37/38
Numerická stabilita
Příklady
► sčítání řady čísel 12 30 + 1 + 1 + 1 + ...je nestabilní
► stabilní alternativa - setřídit a začít od nejmenšího
► úvodní příklad En = 1 - nEn-\
► pro dostatečně velké N položíme Em = 0
► počítáme
F 1
En-1 = -
n
► v každém kroku je chyba redukována faktorem 1/n
Shrnutí
► matematické modely reality
► použití idealizovaného aparátu reálných čísel
► konečná reprezentace čísel je zdrojem problémů
► různé projevy pro různé operace
► míra závadnosti formalizována pojmem numerické stability
► cílem je najít numericky stabilní metodu
► rychlost kritické části výpočtu
► změření nemusí být triviální
► příště se pustíme do konkrétních problémů