01: Úvod do vědeckého počítání

Tento kurz si klade za cíl seznámit studenty s programovacím jazykem Julia.

Julia je moderní programovací jazyk vyvíjený od roku 2012 na MIT. Jde o vysoko-úrovňový dynamický jazyk snažící se dosáhnout vysokého výkonu pomocí JIT kompilace využívající kompilátoru LLVM a možnosti typové anotace.

Komunita kolem Julia se nejvíce zaměřuje na numerické a vědecké výpočty. Z toho pramení i různé vlastnosti návrhu Julia, jako například možnost použití matematických unicode symbolů přímo v kódu, podpora reproducibility výpočtů pomocí správy závislostí v rámci projektu, elegantní podpora práce s maticovými objekty aj. V současnosti ale Julia přetéká i do dalších "více IT" odvětví, datových věd a strojového učení.

Julia je distribuována pod MIT licencí a její zdrojový kód je volně přístupný na Githubu. Instalací Julia a způsoby jak s ní pracovat se budeme zabývat v příští lekci.

O motivaci autorů Julia se lze podrobněji dočíst v poutavém článku [1].

Vědecké výpočty (Scientific Computing)

Vzhledem k slabé expozici studentů FITu fyzikálním problémům se nyní pokusíme nastínit, co přesně vlastně máme na mysli pod "vědeckými výpočty" a jak tato problematika souvisí s historickým vývojem informatiky.

Vznik prvních elektronických počítačů (z úvah vynecháme ty mechanické a živé) spadá do první poloviny (spíše jejímu konci) dvacátého století. Pokud si dohledáte problémy, které se pomocí počítačů lidé snažili vyřešit, tak zjistíte, že velmi často šlo o řešení výpočetních fyzikálních problémů. Řada těchto problémů nemá analytické řešení, lze nalézt jenom přibližná řešení a to za cenu velkého "chroustání čísel" (number crunching).

Za všechny tyto problémy uveďme alespoň následující:

ENIAC (jeden z prvních elektronických počítačů) byl sestrojen (mimo jiné) pro výpočet balistických křivek pro dělostřelectvo.
Matematik John von Neumann, kterého znáte pro po něm pojmenovanou počítačovou architekturu, byl ústředním členem týmu z laboratoře v Los Alamos pracujícím na vodíkové bombě. První počítače využíval právě pro simulace jaderných reakcí.
Alan Turing inicioval vývoj elektronických počítačů při lámání německých šifer (Enigma) během druhé světové války.
První počítače byly hned od začátku používány k numerickému simulování počasí a zejména k jeho předpovídání.

Je až fascinující, co vše mají tyto problémy (a snaha je řešit) v dnešní době za následky a využití. Průměrný dnešní mobilní telefon je mnohonásobně výkonnější než první počítače a výše uvedení velikáni by žasli nad tím, co nosíme bezmyšlenkovitě po kapsách.

Pod vědeckými výpočty tedy máme zhruba na mysli výpočty motivované přírodními jevy, při kterých musí stroj provádět velké množství (ne nutně jenom) aritmetických operací za účelem řešení nějakého (ne nutně jenom) praktického problému.

Programovací jazyky

Programovat první počítače nebylo vůbec jednoduché, používaly se tzv. děrné štítky. Program tedy vlastně byl hromada seřazených štítků, které se samozřejmě nesměly pomíchat.

První programovací jazyky v dnešním smyslu se objevily až v druhé polovině dvacátého století. Jak jsme zmínili dříve, vědecké výpočty typicky vyžadují provádění velkého množství operací. Z toho důvodu je klíčová efektivita a rychlost programu. Proto se pro tyto výpočty historicky používá nejčastěji FORTRAN nebo C/C++.

O C/C++ jste jistě slyšeli minimálně během svého dřívějšího studia. Počátky C můžeme zařadit do sedmdesátých let dvacátého století. Mimo tuto zmínku se mu nyní podrobněji věnovat nebudeme.

FORTRAN

S FORTRANem už jste pravděpodobně do styku nepřišli. Kupodivu jde o jeden z úplně prvních programovacích jazyků, vyvinulo ho IBM už v padesátých letech dvacátého století. Hned od začátku byl zamýšlen jako programovací jazyk pro vědecké a inženýrské výpočty. Název FORTRAN odkazuje na slovní spojení FORmula TRANslation system.

Poslední revize FORTRANu je z roku 2018! Asi se ptáte, jak je možné, že je tu dodnes? Je v něm totiž napsána řada kódu na kterém stojí další vrstvy a technologie. Jako konkrétní ukázku tohoto efektu zmiňme knihovnu BLAS (Basic Linear Algebra Subprograms). Zdrojový kód je samozřejmě napsaný ve FORTRANu. Knihovnu BLAS objevíte například ve vnitřnostech Numpy (Python knihovna pro maticové a numerické výpočty), nebo i v Mathematica, MATLABu, aj.

Jako malou ochutnávku k zamyšlení laskavému čtenáři nabídneme fragmenty kódu z zlaqr4.f (problematika výpočtu vlastních čísel matice; Course pages bohužel neumí správně zvýraznit syntaxi FORTRANu):

*        ==== Hope for the best. ====
*
         INFO = 0
*
*        ==== Set up job flags for ILAENV. ====
*
         IF( WANTT ) THEN
            JBCMPZ( 1: 1 ) = 'S'
         ELSE
            JBCMPZ( 1: 1 ) = 'E'
         END IF
         IF( WANTZ ) THEN
            JBCMPZ( 2: 2 ) = 'V'
         ELSE
            JBCMPZ( 2: 2 ) = 'N'
         END IF
*
*        ==== NWR = recommended deflation window size.  At this
*        .    point,  N .GT. NTINY = 11, so there is enough
*        .    subdiagonal workspace for NWR.GE.2 as required.
*        .    (In fact, there is enough subdiagonal space for
*        .    NWR.GE.3.) ====
*
         NWR = ILAENV( 13, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
         NWR = MAX( 2, NWR )
         NWR = MIN( IHI-ILO+1, ( N-1 ) / 3, NWR )
*
*        ==== NSR = recommended number of simultaneous shifts.
*        .    At this point N .GT. NTINY = 11, so there is at
*        .    enough subdiagonal workspace for NSR to be even
*        .    and greater than or equal to two as required. ====
*
         NSR = ILAENV( 15, 'ZLAQR4', JBCMPZ, N, ILO, IHI, LWORK )
         NSR = MIN( NSR, ( N+6 ) / 9, IHI-ILO )
         NSR = MAX( 2, NSR-MOD( NSR, 2 ) )

nebo

*           ==== Select deflation window size:
*           .    Typical Case:
*           .      If possible and advisable, nibble the entire
*           .      active block.  If not, use size MIN(NWR,NWMAX)
*           .      or MIN(NWR+1,NWMAX) depending upon which has
*           .      the smaller corresponding subdiagonal entry
*           .      (a heuristic).
*           .
*           .    Exceptional Case:
*           .      If there have been no deflations in KEXNW or
*           .      more iterations, then vary the deflation window
*           .      size.   At first, because, larger windows are,
*           .      in general, more powerful than smaller ones,
*           .      rapidly increase the window to the maximum possible.
*           .      Then, gradually reduce the window size. ====
*
            NH = KBOT - KTOP + 1
            NWUPBD = MIN( NH, NWMAX )
            IF( NDFL.LT.KEXNW ) THEN
               NW = MIN( NWUPBD, NWR )
            ELSE
               NW = MIN( NWUPBD, 2*NW )
            END IF
            IF( NW.LT.NWMAX ) THEN
               IF( NW.GE.NH-1 ) THEN
                  NW = NH
               ELSE
                  KWTOP = KBOT - NW + 1
                  IF( CABS1( H( KWTOP, KWTOP-1 ) ).GT.
     $                CABS1( H( KWTOP-1, KWTOP-2 ) ) )NW = NW + 1
               END IF
            END IF
            IF( NDFL.LT.KEXNW ) THEN
               NDEC = -1
            ELSE IF( NDEC.GE.0 .OR. NW.GE.NWUPBD ) THEN
               NDEC = NDEC + 1
               IF( NW-NDEC.LT.2 )
     $            NDEC = 0
               NW = NW - NDEC
            END IF

Moderní nástroje

Postupem času se vedle "nízko-úrovňových" jazyků (jako jsou C a FORTRAN) pro výpočty ve vědecké komunitě čím dál tím více začalo využívat jazyků a nástrojů, u kterých uživatel pracuje na vyšší úrovni. Například Mathematica, MATLAB a jeho open-source alternativa Octave, R, ekosystém balíčků pro Python, aj. Objevily se i možnosti zvýšení výkonu pomocí paralelizace výpočtů. Tyto moderní nástroje většinou využívající dynamické typování, nabízí uživateli větší komfort a pohodlnější pracovní prostředí. Uživatel za to však často musí zaplatit cenu ve formě nižšího výpočetního výkonu.

Tím často vzniká tzv. problém dvou jazyků, kdy prvotní experimenty a prototypy jsou vytvořeny třeba v Mathematica, ale k jejich finální realizaci je nutné řešení přepsat do C nebo FORTRANu (nebo jiného nízkoúrovňového jazyka).

Následující obrázek ukazuje výsledky několika benchmarků jednoduchých výpočetních problémů řešených v různých prostředích relativně vzhledem k C.

Porovnání rychlosti různých jazyků a prostředí. — Obrázek 1. Porovnání rychlosti implementace různych benchmarků relativně vzhledem k implementaci v C [Převzato z Julia homepage].

Řádové rozdíly ve výsledcích jsou zásadní! Pokud provádíte předpověď počasí na zítřejší den, tak není akceptovatelné dostat výsledek za týden! Skutečně se bavíme i o výpočtech, které mohou běžet hodiny, dny a týdny na i velmi výkonných počítačových systémech. Zde každé zlepšení efektivity výpočtu stojí za úsilí.

Proč Julia?

Cíle Julia nejsou malé. Vlastně jde o pokus aplikovat nové prvky (vlastnosti moderních programovacích jazyků, moderní infrastruktura,…) i na oblast vědeckých výpočtů, kde doposud vládne C a FORTRAN. To vše ovšem bez podstatné ztráty výkonu!

Je to utopie? V tento okamžik ještě těžko říci. Uvidíme, jaký budete mít názor na konci tohoto semestru.

Reference

Julia: A Fresh Approach to Numerical Computing. Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah. (2017) SIAM Review, 59: 65–98. doi: 10.1137/141000671.