Fakulta elektrotechnická Katedra matematiky
Dokumentace k semestrální práci
Implementace numerických metod v jazyce C a Python 2013/14 Michal Horáček a Petr Zemek
Vyučující: Mgr. Zbyněk Vastl Předmět: Matematika pro FEL 4
Obsah 1 Zadání semestrální práce 1.1 Název . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Zadání . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 3 3
2 Analýza úlohy 2.1 Ošetření uživatelského vstupu . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Zjištění vlastností matic . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Kontrola špatně podmíněných matic . . . . . . . . . . . . . . . . . . . .
4 4 4 4
3 Matematický popis jednotlivých úloh 3.1 Gausova eliminace . . . . . . . . 3.2 Jacobiho metoda . . . . . . . . . 3.3 Metoda Gausse-Seidela . . . . . . 3.4 Metoda sdružených gradientů . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
5 5 5 6 6
4 Popis implementace v C 4.1 Ošetření vstupu . . . . . . . . 4.2 Datové typy a správa paměti 4.3 Přenositelnost . . . . . . . . . 4.4 Debugování kódu . . . . . . . 4.5 Přehled návratových kódů . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
7 7 7 8 8 8
. . . . .
. . . . .
5 Popis implementace v Pythonu 6 Uživatelská příručka programu 6.1 Překlad a spuštění . . . . 6.2 Chybové hlášky . . . . . . 6.3 Nápověda . . . . . . . . . 6.4 Formát matice soustavy .
9 v . . . .
jazyce . . . . . . . . . . . . . . . .
C . . . . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
10 10 10 11 11
7 Závěr 7.1 Výkonnost programu . . . . . . . . . . 7.1.1 Program v jazyce C . . . . . . 7.1.2 Program v jazyce Python . . . 7.2 Možný problém na určitých systémech 7.3 Vysvětlení nestandardních konstrukcí . 7.4 Kódování zdrojových souborů . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
12 12 12 13 13 14 14
2
1 Zadání semestrální práce 1.1 Název Implementace numerických metod v jazyce C a Python
1.2 Zadání Naprogramujte v jazyce C a Python přenositelnou konzolovou aplikaci, která jako vstup načte z parametru na příkazové řádce jméno souboru s maticí soustavy a ostatní potřebné údaje pro spuštění výpočtu. Implementujte Gausovu eliminační metodu, metodu sdružených gradientů, Jacobbiho metodu a Gauss-Seidelovu metodu řešení soustav linárních algebraických rovnic. Implementaci v jazyce C řeší Petr Zemek a v jazyce Python Michal Horáček.
3
2 Analýza úlohy 2.1 Ošetření uživatelského vstupu Představme si, co všechno může BFU uživatel zkoušet zadávat do vstupu programu. Je vhodné ošetřit případy, kdy je očekáváno číslo. Pokud uživatel zadá nečíselný znak, měla by se použít defaultní hodnota. Matice je vhodné načítat spíše ze souboru než z výchozího vstupu pro jejich rozsáhlost. Program by měl umět parsovat soubor s maticí soustavy. Následně by měl umět ověřit, jestli matice v souboru je opravdu čtvercová a nechybí-li některý člen.
2.2 Zjištění vlastností matic Pro celistvost programu by bylo vhodné zjišťovat vlastnosti matic a na základě získaných informací se rozhodnout, jaký algoritmus pro řešení zvolit.
2.3 Kontrola špatně podmíněných matic Před spuštěním algoritmu Jacobiho metody nebo Metody Gausse-Seidela je vhodné ověřit, že se v matici soustavy na hlavní diagonále nenachází nula, protože se v těchto metodách vektor neznámých dělí prvkem na hlavní diagonále.
4
3 Matematický popis jednotlivých úloh 3.1 Gausova eliminace Matici soustavy je třeba převést do stupňovitého tvaru na horní trojúhelníkovou matici. Pod hlavní diagonálou musí zbýt nuly. Čili na konci eliminace musí zbýt na posledním řádku a sloupci matice jedno číslo, které po vynásobení neznámou z posledního řádku a sloupce se rovná poslednímu řádku vektoru pravé strany. Pokud při eliminaci vyjde nulový řádek, soustava má nekonečně mnoho řešení. Obecný algoritmus: for k = 1 ... m: #Najít sloupec, kde je největší číslo (pivot): i_max := argmax (i = k ... m, abs(A[i, k])) if A[i_max, k] = 0 error "Matice je singulární!" swap rows(k, i_max) # pro každý prvek pod pivotem: for i = k + 1 ... m: # pro všechny zbývající prvky v daném sloupci for j = k + 1 ... n: A[i, j] := A[i, j] - A[k, j] * (A[i, k] / A[k, k]) endfor # Naplnit spodní trojúhelníkovou matici nulami: A[i, k] := 0 endfor endfor
3.2 Jacobiho metoda Jacobiho metoda je iterační metoda, která konverguje pouze pro ryze řádkově nebo sloupcově diagonálně dominantní matice.
x𝑘+1 = 𝑖
b𝑖 A𝑖𝑖
−
𝑛 ∑︀ 𝑗=𝑖+1
A𝑖𝑗 𝑘 x , A𝑖𝑖 𝑖
𝑘 ≥ 0;
𝑖 = 1, ..., 𝑛
5
3.3 Metoda Gausse-Seidela Je téměř shodná s Jacobiho metodou. Liší se v tom, že v každé iteraci se od vektoru neznámých odečítá aproximace z předchozího kroku. Tím se sníží počet potřebných iterací k dosažení stejné přesnosti jako při použití Jacobiho metody.
x𝑘+1 = 𝑖
b𝑖 A𝑖𝑖
−
𝑖−1 ∑︀ 𝑗=1
A𝑖𝑗 𝑘+1 x A𝑖𝑖 𝑖
−
𝑛 ∑︀ 𝑗=𝑖+1
A𝑖𝑗 𝑘 x , A𝑖𝑖 𝑖
𝑘 ≥ 0;
3.4 Metoda sdružených gradientů Obecný matematický zápis:
𝜀> 0, x1 , 𝑘 = 1 b1 = Ax1 − b,
v1 = r1 ,
𝛼1 =
||b1 ||2 , Ar1 , r1
while || x𝑘+1 − x𝑘 || >
𝜀 k=k+1 r𝑘 = r𝑘−1 + 𝛼𝑘−1 Av𝑘−1 𝑘
𝑘−1
r ,Av 𝛽𝑘−1 = − Av 𝑘−1 𝑘−1 ,v
v𝑘 = r𝑘 + 𝛽𝑘−1 v𝑘−1 r𝑘 ,v𝑘−1 Av𝑘 ,v𝑘 𝑘+1 𝑘
𝛼= − x
= x + 𝛼𝑘 v𝑘
end while
6
x2 = x1 + 𝛼1 v1
𝑖 = 1, ..., 𝑛
4 Popis implementace v C 4.1 Ošetření vstupu Nejdříve se ověřují parametry příkazové řádky - jestli uživatel zadal cestu k matici a číslo metody řešení (soubor rozcestnik.c; ř. 29-68). Potom se ověří, jestli existuje soubor s maticí a uživatel, pod kterým běží tento program, má oprávnění ke čtení souboru (rozcestnik.c; ř. 81-85). V případě, že uživatelský vstup prošel těmito kritérii, zjistí se počet sloupců matice(rozcestnik.c; ř. 102-160) a provede se alokace množství paměti, které odpovídá druhé mocnině zjištěných sloupců matice bez pravé strany plus jednoho sloupce navíc (pravá strana). Nyní je vytvořena matice, která se musí naplnit daty ze souboru. To řeší funkce nacist_matici() v souboru jednoduche_fce/nacist_soubor.c. Zde se provádí rozřazování prvků matice ze souboru do matice v paměti a zároveň se provádí jednoduché testování na přítomnost znaku "|" oddělující matici soustavy od pravé strany (jednoduche_fce/nacist_soubor.c; ř. 17-62). Nakonec se porovná počet řádků a sloupců. Pokud počet nesouhlasí, matice není čtvercová a program skončí s návratovým kódem 1.
4.2 Datové typy a správa paměti Nejpodstatnějšími proměnnými v programu jsou matice soustavy a ostatní vektory. Program byl vyvíjen na 32-bitové architektuře a proto pro tyto proměnné byl zvolen datový typ float. Práce s tímto datovým typem na této architektuře je optimální z hlediska rychlosti běhu programu. Pokud bychom chtěli na stejné architektuře např. zdvojnásobit přesnost výpočtu, můžeme zvolit datový typ matic a vektorů s dvojnásobou přesností (double float) a nebo dvojnásobně zmenšit 𝜀. V prvním případě docílíme toho, že se jedno desetinné číslo bude kopírovat dva procesorové cykly dlouho. Matematické operace mohou trvat až několikanásobně více procesorových cyklů déle než při použití jednoduché přesnosti. V druhém případě sice vzroste režije algoritmu, která se skládá pouze z inkrementace proměnných typu unsigned int a udržovacích podmínek, ale hlavně vzroste počet iterací, které jsou závislé na použité metodě. Efektivnější je tedy buď změnit 𝜀, a nebo použít procesor s širší sběrnicí a přepsat konstantu datového typu v souboru headers/konstanty.h. Paměť je z většiny přidělována staticky. Ta, která se alokuje dynamicky by se mohla alokovat staticky, protože v průběhu běhu programu ji není kde dealokovat. Paměť se alokuje na začátku funkce main() (ř. 131-160) hned po zjištění rozměrů matice soustavy. Dealokuje se na konci programu - ať už skončí chybou nebo ne.
7
4.3 Přenositelnost Program není závislý na softwarové ani hardwarové platformě a neobsahuje ani absolutní cesty k souborům. Jako parametr lze předávat relativní i absolutní cesty. Pro jednodušší úpravu zdrojových kódů jsou použity konstanty. Při kompilaci programu pro HW platformu s šířkou sběrnice vyšší než 32 bitů je vhodné ověřit, jestli float má skutečně tolik bitů jako šířka sběrnice.
4.4 Debugování kódu ...se provádí prostým odkomentováním řádky #define _debug v souboru headers/konstanty.h.
4.5 Přehled návratových kódů Návratový kód 0 1 2 3 4 5 6 7 10 11 12 20 21 22 23 24
Význam Program skončil úspěšně Uživatel nezadal číslo metody řešení Uživatel nezadal cestu k souboru s maticí Soubor s maticí soustavy neexistuje, nebo uživatel nemá oprávnění ke čtení Operační systém odmítl přidělit paměť pro matici soustavy (není žádná volná) Uživatel zvolil chybné číslo metody řešení Chyba syntaxe v matici soustavy: očekáván znak "|" Uživatel zadal neznámý parametr příkazu Chyba v programu; funkce gaus_seidel(), soubor vypocetni_fcegaus-sidel.c Chyba v programu; funkce jacobi(), soubor vypocetni_fcejacobi.c Chyba v programu; funkce sdruzene_gradienty(), soubor vypocetni_fcegrad.c Výpočet metodou Gausse-Seidela ukončen pro dosažení maximálního počtu iterací Výpočet metodou Sdružených gradientů ukončen pro dosažení maximálního počtu iterací Výpočet Jacobiho metodou ukončen pro dosažení maximálního počtu iterací Při pokusu o spuštění algoritmu Jacobiho metody byla na diagonále matice soustavy zjištěna nula Při pokusu o spuštění algoritmu Gausse-Seidela byla na diagonále matice soustavy zjištěna nula
8
5 Popis implementace v Pythonu Hlavní výhodou je jednoduchost implementace. Jazyk Python je přímo navržen pro snadné a efektivní psaní aplikací, tudíž programátor neřeší nízkoúrovňové operace a může se plně soustředit na hlavní účel aplikace. Jelikož je to však interpretovaný jazyk, je nutné počítat s nižší rychlostí výpočtů, avšak v případě Pythonu není ztráta výkonu tak velká, neboť výkonově kritické knihovny jsou napsány v jazyce C. Nabízí se otázka, zda v Pythonu implementovat metody pro řešení soustav lineárních rovnic. Z vlastností jazyka plyne, že by nebylo vhodné jej používat pro velké soustavy rovnic, nebo opakovaně spouštět výpočet. V druhém případě by nastal problém, kdy se spotřebuje určitý čas pro spuštění interpretru Pythonu, a teprve potom se spustí výpočet samotný. Při použití všech čtyř metod na soustavu 2x2 trval běh výpočtu (vč. spuštění interpretru) 41 ms. Spuštění samotného interpretru trvalo 30 ms, tudíž výpočet zabral 11 ms, a to pro všechny čtyři metody. To znamená, že pro tak malou soustavu bylo spotřebováno přibližně 73% času jenom pro rozběhnutí výpočtu. Pokud by bylo potřeba počítat malé soustavy opakovaně, bylo by třeba interpretr spustit jednou a poté v něm opakovaně spouštět výpočty. Nepřekvapí ani paměťová náročnost. Na konci výpočtu proces využívá zhruba 3.5 MB RAM. Naštěstí v Pythonu existuje Garbage Collector, který nepotřebné proměnné dealokuje z paměti, takže i když programátor nedbá na paměťovou náročnost, neroste spotřeba paměti s každou proměnnou. Další výhodou implementace v Pythonu je skutečnost, že lze data z výpočtů snadno dále zpracovávat, například s pomocí knihovny matplotlib lze vytvářet grafy jako v prostředí Matlab. Interpretr Pythonu lze také snadno integrovat do jiných aplikací, kde jsou výpočty potřeba (takovou “konzoli” lze vidět například v aplikaci pro 3D modelování, Blender, kde o výpočty není nouze). Výhodou jsou také snadno dostupné funkce knihovny pro výstup, ať už v textové podobě, HTML, LATEXu nebo PDF.
9
6 Uživatelská příručka programu v jazyce C 6.1 Překlad a spuštění Pro překlad bez zásahu do zdrojových kódů je nutné mít dostupný překladač gcc. Pro samotný překlad stačí zadat příkaz make v kořenovém adresáři zdrojových kódů (nachází se tam soubor Makefile. Program má jednoduché textové rozhraní. Spouští se z příkazové řádky; standardní rozměry konsole 80x25 postačují. Pro spuštění je třeba zadat cestu k matici a metodu řešení. Příklad:
./slar -s matrix3.txt -m 0 Tím dojde k načtení souboru matrix3.txt a použití metody 0 - GEM. Parametr -m může nabývat těchto hodnot: Hodnota 0 1 2 3
Použitá metoda Gausova eliminační metoda Metoda Gausse-Seidela Metoda Sdružených gradientů Jacobbiho metoda
V závisloti na použité metodě můžete být požádání o další uživatelský vstup - o hodnotu 𝜀 a o počáteční vektor neznámých. Více informací nebude vyžadováno. Pokud zadáde jiný znak než číslo, použije se výchozí hodnota. Výchozí hodnota 𝜀 je 0.01 a vektor neznámých je defaultně 0.
6.2 Chybové hlášky Program při svém ukončení předává operačnímu systému návratový kód. Podle tohoto kódu lze usoudit, zda program skončil správně nebo chybně. Při svém ukončování také vypíše, proč končí. Např. když se stane, že ve vašem PC už nezbývá dostatek volné paměti (používáte příliš novou verzi Windows), vypíše se hláška: Chyba #4: Došla paměť!
10
6.3 Nápověda Pro výpis nápovědy zadejte: ./slar -h
6.4 Formát matice soustavy Do načítaného souboru se zapisuje rozšířená matice soustavy. Matice A musí být čtvercová. Jako oddělovač pravé strany se používá znak "|". Jako oddělovač sloupců se používá mezera. Před a za oddělovačem pravé strany musí být mezera. Příklad obsahu souboru s rozšířenou maticí soustavy: 1 3 4 | 1 8 4 3 | 8 9 2 4 | 9
11
7 Závěr 7.1 Výkonnost programu Výkonnost programu v C a programu v Pythonu se ověřovala na stroji IBM T34 s procesorem Intel Pentium-M s taktovací frekvencí 18644 MHz. L1 cache procesoru je 64kiB, L2 cache 2MiB. Šířka pamětí, procesoru a FSB sběrnice je 32 bitů. Operační systém je linuxová distribuce Gentoo. Testovaný program byl spuštěn s prioritou nice -10, aby se omezily vlivy ostatních aplikací. Operační systém běžel pouze v příkazovém režimu, aby se redukoval počet spuštěných procesů. Program psaný v jazyce C byl testován na rychlost bez debugovacích funkcí a s defaultní hodnotou 𝜀, maximálního počtu iterací a počáteční hodnotou vektoru neznámých. Toho se docílilo zakomentováním části kódu, která obstarává výchozí vstup. Matice použitá k měření má rozměry 200 x 200 a je přiložená v souboru matice_velka.txt v kořenovém adresáři zdrojových kódů.
7.1.1 Program v jazyce C GEM Čas běhu: 243 ms
Jacobiho metoda Čas běhu: 87 ms Dosažená přesnost: 0.008262 Počet iterací: 25
Medoda sdružených gradientů Čas běhu: 122 ms Dosažená přesnost: 0.009407 Počet iterací: 9
12
Metoda Gausse-Seidela Čas běhu: 60 ms Dosažená přesnost: 0.003813 Počet iterací: 17
7.1.2 Program v jazyce Python GEM Čas běhu: 6314.017 ms
Jacobiho metoda Čas běhu: 1170.088 ms Počet iterací: 35
Medoda sdružených gradientů Čas běhu: 1179.550 ms Počet iterací: 7
Metoda Gausse-Seidela Čas běhu: 248.547 ms Počet iterací: 9
7.2 Možný problém na určitých systémech Je zde jedna závada, na kterou bych chtěl upozornit - jedná se o uvolňování paměti. Objevuje se pod 32-bitovým Gentoo (i686), pod 32-bitovým Archlinuxem (arm-v6j) i na 32-bitové ditribuci Fedora 18 (i686) (kernel 3.11.10-100.fc18.i686 s glibc 2.1634) Při použití matic menších než 7x7 se chyba neobjevuje vůbec. Na 64-bitových architekturách se paměť dealokuje v pořádku. Nástroj Valgrid vypíše, že byla veškerá paměť uvolněna, i když dealokace skončí s chybou. Výpis chyby: *** glibc detected *** ./slar: free(): invalid next size (fast): 0x0804c200 *** ======= Backtrace: ========= /lib/libc.so.6(+0x79930)[0xb7753930] ./slar[0x80487b4] ./slar[0x8048bcd]
13
/lib/libc.so.6(__libc_start_main+0xf1)[0xb76f6c31] ./slar[0x80486bd] ======= Memory map: ======== 08048000-0804b000 r-xp 00000000 08:06 170979 /home/jashin/Dokumenty/Vyvoj/Rozpracovane/M4E/slar ... ...
7.3 Vysvětlení nestandardních konstrukcí V programu jsem využil nekonečných cyklů a příkazů skoku. Nekonečné cykly jsou nekonečné proto, že se tím zjedoduší programování a zpětná analýza kódu. Ukončovací podmínky jsou někde uvnitř cyklů a tudíž nemůže být použit např. cyklus for, který má udržovací podmínku v parametru. V kódu jsem příliš nepoužíval standardní ANSI-C komentáře - při vývoji jsem někdy zakomentovával velké kusy kódu a zakomentovaný komentář ( /* /*komentář */ */ ) se překladačem vyhodnotí jako chyba. Proto používám komentáře z C++.
7.4 Kódování zdrojových souborů Všechny zdrojové kódy programu i dokumentace jsou uloženy v kódování UTF-8 s UNIXovým ukončením řádků.
14
