C Praktické základy numerických výpočtů
Smyslem této kapitoly není podávat systematický popis základních metod numerické matematiky, v tomto směru odkazuji zájemce například na skriptum Humlíček (2009) nebo odpovídající učebnice (například Přikryl, 1985; Vitásek, 1987; Čermák & Hlavička, 2006), ale pouze stručně a názorně ilustrovat některé principy a možné postupy při praktickém numerickém modelování nejčastěji se v praxi vyskytujících (a výše popsaných) analytických okruhů. Vybrané příklady jednoduchého modelování jsou také doprovázeny ukázkami programovacích skriptů pro daný konkrétní problém, případně obrázky a grafy výsledných modelů.
Rozsáhlé využívání numerické matematiky ve většině přírodovědných a technických disciplín přineslo také tvorbu celé řady hotových knihoven, rozepsaných do hlavních programových jazyků; jejich přehledné úložiště se nachází například na stránkách GAMS (Guide to Available Mathematical Software). Některé z nich jsou komerční (a poměrně komplexní), například NAG (Numerical Algorithmus Group) nebo IMSL (International Mathematics and Statistics Library), jiné jsou volně dostupné a bývají zpravidla zaměřené na specifickou oblast, například FFTPACK – (rychlá) Fourierova transformace, LAPACK (viz odstavec C.1) – lineární algebra, MINPACK – nelineární rovnice, atd. Jejich úplná nebo i částečná implementace do vlastních algoritmů může výrazně urychlit a usnadnit jejich tvorbu i kvalitu.
C.1 Numerické metody lineární algebry
V tomto oddíle nebudeme probírat jednotlivé metody numerických řešení lineární algebry, k základům tohoto tématu existuje rozsáhlá literatura (např. Humlíček (2009)), popisující jak vlastní numerické rovnice tak jejich stabilitu a podmíněnost (tj. zejména stanovení přesnosti numerických maticových algoritmů), odhady chyb, atd. V současnosti existuje řada hotových balíčků (procedur), sestávajících z jednotlivých podprogramů (knihoven), určených pro řešení dílčích nebo i kombinovaných algebraických úloh (například řešení soustav lineárních rovnic, řešení tzv. tridiagonálních matic (tj. matic s nenulovými prvky pouze na hlavní a obou sousedních diagonálách), hledání determinantů, inverzních matic, vlastních hodnot a vlastních vektorů, atd. Jedním z nejvýkonnějších takových balíčků, který zde podrobněji představíme, je programový balíček LAPACK (Linear Algebra PACKage, Anderson et al. (1999)), který se vyvinul ze starších balíčků EISPACK a LINPACK a je určen pro fortran 77, fortran 90, existují také C++ verze. Existují rozšířené verze tohoto balíčku nebo další knihovny na něm postavené, se zabudovanými podprogramy pro paralelizaci na výkonných počítačových clusterech (viz odstavec C.5), například ScaLAPACK, MAGMA, MORSE, CHAMELEON, atd.
Programový balíček LAPACK je volně dostupná softwarová knihovna, její instalaci provedeme buď ze softwarového centra používané
systémové distribuce, nebo z adresy http://www.netlib.org/lapack. Při překladu (kompilaci) používaného programového souboru zadáme odkaz na
LAPACK, například: gfortran název_souboru.f95 -llapack. Popis jednotlivých podprogramů a jejich využití
(např. knihovna DGBSV pro řešení
soustav reálných lineárních rovnic o libovolném počtu proměnných nebo DGTSV, vhodná pro řešení tridiagonálních matic, atd.) jsou dostupné v uživatelských příručkách, například v Anderson et al. (1999).
| N | (vstup, INTEGER) = počet rovnic = řád čtvercové matice A, N ≥ 0. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| KL | (vstup, INTEGER) = počet spodních diagonál matice A, KL ≥ 0. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| KU | (vstup, INTEGER) = počet horních diagonál matice A, KU ≥ 0. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| NRHS | (vstup, INTEGER) = počet sloupců pravé strany (tj. matice B), NRHS ≥ 0. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| AB | (vstup/výstup, DOUBLE PRECISION) = pole dimenze (LDAB,N). Na vstupu: matice A v pásovém uložení, v řádcích od KL+1 do 2*KL+KU+1; řádky 1 až KL pole nemusí být vypsány. j-tý sloupec pole A je uložen jako j-tý sloupec pole AB následovně: AB(KL+KU+1+i-j, j) = A(i,j) pro max(1, j-KU) ≤ i ≤ min(N, j+KL). Na výstupu: detaily faktorizace – matice U je uložena jako horní trojúhelníková pásová matice s KL+KU horními diagonálami v řádcích od 1 do KL+KU+1, multiplikátory M, použité během faktorizace jsou uchovány v řádcích od KL+KU+2 do 2*KL+KU+1 (viz schéma níže). |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| LDAB | (vstup, INTEGER) = určující dimenze pole AB. LDAB ≥ 2*KL+KU+1. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| IPIV | (výstup, INTEGER) = pole dimenze (N), indexy pivotů, které definují permutační matici; i-tý řádek matice byl zaměněn za řádek IPIV(i). | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| B | (vstup/výstup, DOUBLE PRECISION) = pole dimenze (LDB,NRHS), na vstupu je N : NRHS matice pravé strany B. Na výstupu, pokud INFO = 0, je N : NRHS řešení matice X. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| LDB | (vstup, INTEGER) = určující dimenze pole B. LDB ≥ max(1,N). | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| INFO | (výstup, INTEGER) = 0: úspěšný výstup, < 0: pokud INFO = -i, pak i-tý argument má nepovolenou hodnotu, > 0: pokud INFO = i, U(i,i) je přesně 0. Faktorizace je ukončena, ale faktor U je přesně singulární, řešení nemohlo být vypočteno. | |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| DALŠÍ DETAILY: | Schéma pásového uložení je ilustrováno na následujícím příkladu, kdy M = N = 6, KL = 2, KU = 1:
Prvky pole označené * nejsou používané ve výpočetním procesu; prvky označené + nemusí být uvedeny na vstupu, ale jsou nutné ve výpočetním procesu pro uložení prvků pole U z důvodu nedostatku místa, vyplývajícího z výměny řádků. |
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| N | (vstup, INTEGER) = počet rovnic = řád čtvercové tridiagonální matice A, N ≥ 0. |
| NRHS | (vstup, INTEGER) = počet sloupců pravé strany (tj. matice B), NRHS ≥ 0. |
| DL | (vstup/výstup, DOUBLE PRECISION) = na vstupu pole prvků spodní (sub) diagonály matice A, dimenze N-1, na výstupu je toto pole přepsáno N-2 prvky druhé horní diagonály horní trojúhelníkové matice U, dané LU faktorizací. |
| D | (vstup/výstup, DOUBLE PRECISION) = pole dimenze N, na vstupu obsahuje diagonální prvky matice A, na výstupu je toto pole přepsáno diagonálními prvky matice U. |
| DU | (vstup/výstup, DOUBLE PRECISION) = pole dimenze N-1, na vstupu obsahuje N-1 prvků horní (super) diagonály matice A, na výstupu je toto pole přepsáno N-1 prvky první horní diagonály horní trojúhelníkové matice U. |
| B | (vstup/výstup, DOUBLE PRECISION) = pole dimenze (LDB,NRHS), na vstupu je N : NRHS matice pravé strany B. Na výstupu, pokud INFO = 0, je N : NRHS řešení matice X. |
| LDB | (vstup, INTEGER) = určující dimenze pole B. LDB ≥ max(1,N). |
| INFO | (výstup, INTEGER) = 0: úspěšný výstup, < 0: pokud INFO = -i, pak i-tý argument má nepovolenou hodnotu, > 0: pokud INFO = i, U(i,i) je přesně 0. Faktorizace je ukončena, ale faktor U je přesně singulární, řešení nemohlo být vypočteno. |
Obdobným způsobem jsou sestaveny i ostatní knihovny. Příklady řešení a programových skriptů s odkazem na LAPACK uvádíme v dalších odstavcích C.2.1, C.2.2, C.2.3, C.2.4, C.4.1, C.4.7, atd.