Fortran: Polní výzbroj

Evoluce

Pole existovala ve Fortranu odedávna. Poskytovala možnost abstrakce od přímé práce s pamětí pomocí ukazatelů. Ještě ve Fortranu 77 byla pole jen datovou strukturou – bylo možné jen číst nebo zapsat jednotlivý prvek, případně předat pole jako argument procedury. V době, kdy se dokončoval Fortran 77, který se měl stát průmyslovým standardem pro výpočetní programy, vznikal na univerzitě v New Mexico (USA) i jiný produkt, který čekala velmi úspěšná budoucnost. Cleve Moler vytvořil Matlab proto, aby svým studentům zpřístupnil fortranské knihovny pro lineární algebru LINPACK a EISPACK, aniž by se museli učit Fortran 77. Za svůj prvotní úspěch tak vděčil Matlab nejen elegantní a interaktivní práci s vektory a maticemi, ale ve stejné míře těmto kvalitním knihovnám (interpret byl původně rovněž ve Fortranu, ale byl poměrně záhy přepsán do C). Matlab ovšem tuto „půjčku“ splatil, neboť není pochyb, že aritmetika polí ve Fortranu 90 je do značné míry „obšlehnutá“ z Matlabu. Prvenství v této oblasti ale nepatří ani Matlabu – už před ním existoval ambiciózní jazyk APL (který se ale příliš nerozšířil zřejmě kvůli kryptické syntaxi a problematické implementaci). Pokud si s Matlabem nebo jeho klony (Octave, Scilab) rozumíte, zřejmě vám bude práce s poli ve Fortranu docela blízká.

Deklarace polí

Jednoduchou deklaraci pole jsme už viděli. Reálný či komplexní vektor o deseti prvcích se deklaruje stylem:

real:: vector1(10)
! nebo
complex,dimension(10):: vector2

Oba způsoby deklarace jsou zcela rovnocenné – vše, co lze uvést do závorek za proměnnou, lze uvést i do atributu dimension a obráceně. První způsob je zřejmě bližší jazyku C a méně „ukecaný“, druhý má naopak blíže k Pascalu a přijde mi trochu logičtější s ohledem na myšlenky předchozích odstavců. Osobně klidně kombinuji oba způsoby podle toho, co se víc hodí. Dlouhá slova jako dimension nepředstavují s „namakanými“ programátorskými editory problém – já například druhou deklaraci výše ve ViMu vytvořím posloupností kláves: co`,dm`10<Ctrl-J>:: vector2.

Jak už také víme, znamená definice dimension(10) rozsah indexů 1 až 10. Můžeme to změnit specifikací dolní mez:horní mez, například na indexování od nuly: dimension(0:9). Povolena jsou i záporná čísla. Fortran umí pracovat i s vícedimen­zionálními po­li:

real:: matrix(m,n),tensor(3,0:m,0:n,2)

Počet dimenzí se ve fortranštině označuje jako „rank“ a může být maximálně 7. Pro každou dimenzi můžeme opět specifikovat libovolný rozsah indexů. Takováto pole jsou „kartézská“, tedy každý index může nabývat nezávisle všech hodnot (některé prvky můžeme samozřejmě ignorovat, ale stejně budou zabírat paměť). Vektor délek rozsahů jednotlivých dimenzí pole se nazývá tvar (shape). Moderní Fortran umožňuje vytvářet i „zubatá“ vícedimenzionální pole, tedy pole polí různých velikostí, ale o tom až jindy. Poznamenejme, že Fortran zplnoprávňuje i pole nulové velikosti – to když je v některé dimenzi prázdný rozsah, např 1:0. Ta sice obvykle nebudete samostatně vytvářet (co taky pak s nimi), ale často vznikají jako mezní případ, a protože jim Fortran rozumí, odpadá tak v mnoha případech jejich zvláštní ošetřování.

Určení mezí a alokace

Podle způsobu určení velikosti pole (a „za oponou“ i alokace, tedy přidělení paměti) lze ve Fortranu rozlišit tři druhy polí: statická,auto­matická a dynamická. (To platí pro deklaraci proměnných, ne formálních argumentů procedur).

• Statická pole mají konstantní meze – tedy jsou určitelná již při překladu. Lze použít samozřejmě jak literály, tak pojmenované konstanty ( parameter) a různé výrazy (ne však všechny, např. množina použitelných funkcí je omezená.) Tato pole mohou být deklarována jako „globální“ proměnné modulů nebo lokální proměnné procedur a také jako konstanty parameter. Obvykle se v paměti umístí do datového segmentu výsledného programu.
• Automatická pole lze použít jen jako proměnné procedur. V deklaraci mezí smíme použít kromě konstant také formální argumenty procedury a dokonce i globální proměnné, ke kterým má procedura přístup. Meze pole se tak vždy určí při vstupu do procedury a po dobu jejího běhu jsou neměnné. Tato pole se typicky alokují na stacku.
• Dynamická pole jsou deklarována s atributem allocatable nebo pointer a místo specifikace mezí mají jen dvojtečku (tzv. deferred shape). Za běhu procedury jim pak lze určit meze a přidělit paměť příkazem allocate a dealokovat je příkazem deallocate. Mohou být jak globální proměnné modulů, tak lokální proměnné procedur. Dealokace lokálních polí s atributem allocatable je navíc automatická, udělá ji za nás překladač. Tato pole, jak asi tušíte, se obvykle umisťují do dynamické paměti – heap.

Příklady:

module pos
implicit none
real:: Q(3,3),t(3)! statická pole
real,allocatable:: pts(:,:),proj(:,:) ! dynamická pole

contains

subroutine do_something(f)
real,intent(in):: f
real:: tpts(3,size(pts,2)) ! automatické pole
real,allocatable:: pp(:) ! dynamicke pole - automaticky dealokované
! ...
end subroutine
end module

Počítání s poli

Základní myšlenka „polní“ aritmetiky je jednoduchá: neuvažovat o poli jako datové struktuře – seznamu hodnot, ale jako o jednoduché proměnné nebo výrazu. Tak běžně pracujeme v matematice s vektory, maticemi a tenzory: jsou to pro nás entity, s nimiž provádíme různé operace. Výhodou tohoto přístupu je zejména jednodušší programování – je často bližší našemu uvažování. Samozřejmě ne vždy – občas naopak potřebujeme uvažovat o poli jako o „úložně“ prvků a je velkým omylem snažit se za každou cenu vyhýbat smyčkám a zpracovávání pole „po prvcích“. Dobrou radou je vždy použít tu konstrukci, která vám momentálně připadá nejčitelnější a nejpřirozenější.

Všechny základní matematické operace a funkce, zmiňované v minulém díle, jsou elementální. To ve Fortranu znamená, že každý operand nebo argument může být jak skalár, tak pole. Operandy, které jsou poli, musejí mít stejný tvar. Výsledkem funkce či operace je pak opět pole, jehož prvky se získají aplikací operace či funkce na vzájemně si odpovídající elementy argumentů a operandů, přičemž skaláry jsou automaticky rozšířeny na pole stejné velikosti. Jsou-li a,b pole stejných rozměrů, x skalár, lze například psát

sqrt(a**2 + b**2)
a * max(b,0.)**x
atan2(a,b+1)

a výsledkem všech těchto výrazů je opět pole stejné velikosti jako a,b.

Přiřazení = funguje i s poli. Lze psát pole = výraz, kde pole je proměnná typu pole (nebo její sekce, viz dále), a výraz se vyhodnotí na skalár nebo pole. V případě pole musí samozřejmě souhlasit rozměry.

Sekce polí

Ve výrazech nebo argumentech procedur nemusíme použít vždy celé pole, je možné specifikovat sekci. To se dělá podobně jako indexování, tedy výběrem prvku, jen místo jednoduchého indexu použijeme takzvaný indexový triplet [dolní mez]:[horní mez][:krok]. Obě meze lze vynechat. V tom případě se použijí vlastní meze pole. Krok lze rovněž vynechat, default je 1. Triplet d:h:k definuje aritmetickou posloupnost indexů d,d+k,...,d+k*((h-d)/k). Mají-li h-d a k opačné znaménko, je posloupnost prázdná. Krok může být i záporný – tedy je možné např. obrátit pole, ale nesmí být nula. V každé dimenzi můžeme specifikovat buď triplet, nebo jednoduchý index. Výsledná sekce je pole, jehož rank je počet specifikovaných tripletů, a rozměry určené délkami posloupností vytvářených triplety. Kartézský součin těchto indexových posloupností spolu s jednoduchými indexy pak určuje prvky původního pole, které tvoří sekci. Na rozdíl od Matlabu je podstatný rozdíl mezi i  a  i:i!.

Kvíz: jaký je výstup následujícího programu? ( shape je funkce, která vrací rozměry pole). Zkuste odpověď nejprve odvodit sami a pak ověřit překladačem.

program kviz
integer:: a(8,-3:5,2:11)
print *,shape( a(2:5,5:-1:-2,3) )
print *,shape( a(1:1,::2,:8:2) )
print *,shape( a(::-1,2,2:) )
end program

Kromě jednoduchého indexu a tripletu lze na stejném místě uvést ještě tzv. vektorový index, což je jednorozměrné pole integerů. Vytváří se opět sekce, jen indexová posloupnost pro příslušnou dimenzi je určena přímo tímto vektorem. Používáme-li pole s vektorovým indexem jako levou stranu přiřazení, nesmí se ve vektoru indexy opakovat, jinak mohou. Pole s vektorovým indexem také nelze použít pro out nebo inout argument procedury (viz dále). Příklady:

integer:: p(m) ! permutace
real :: a(m,n),b(m)
! ...
a = a(p,:) ! permutovat radky matice
b(p) = b ! inverzní permutace

Odbočka: Nulaři a jedničkáři

Prvotní specifikace Fortranu umožňovaly deklarovat jen pole indexované od jedničky, jak bylo zvykem od nepaměti. Byla-li p adresa prvního prvku, pak se adresa i-tého prvku vypočítala jako p-1+i. Při výpočtu adresy tak byla potřeba jedna instrukce navíc (dekrementace). To bylo samozřejmě trnem v oku všem, kteří se snažili o to, aby byly výpočty na počítačích efektivní. Záhy se nalezla dvě možná řešení:

1. Indexovat vždy od nuly.
2. Vycházet z adresy (neexistujícího) prvku s indexem nula. To umožňuje efektivní indexování s libovolnou dolní mezí.

Indexování od nuly bylo popularizováno zejména jazykem C, který pole neabstrahuje jako Fortran, ale odkrývá přímo ukazatelovou aritmetiku, takže indexování od nuly je v něm přirozené. Postupně bylo vývojáři a softwarovými inženýry přijato jako standard a dostalo se i do dalších jazyků: C++, Java, Python, Ruby, Perl, Scheme, PHP, JavaScript, Oberon aj. S indexováním od nuly se typicky pojí zadávání rozsahů jako „upper exclusive“, tedy nezahrnujících horní mez – např. C++ STL, Python, Ruby. Výjimkami jsou např. Perl, IDL a Basic (Visual Basic, FreeBasic), používající „upper inclusive“.

Na „druhé straně barikády“ jsou jazyky, nástroje a prostředí, které se drží tradičního indexování od jedničky a „upper inclusive“ rozsahů. Jsou to zejména „méně programátorské“ jazyky a prostředí, např. Fortran, Matlab/Octave/Sci­lab, R, Excel, Mathematica, Maple a další, které nejsou zrovna určeny ani tak pro vývojáře, ale spíš pro vědce, inženýry, ekonomy, účetní apod. Některé jazyky umožňují používat libovolnou dolní mez polí (Fortran, Pascal, FreeBasic, Ada) jakožto rozumný kompromis. To je ovšem u jazyků, které nedeklarují proměnné, trochu problém.

Spor o „správný“ způsob indexování a zapisování rozsahů je zřejmě jedna z vleklých „flamewars“ moderní doby. Indexování od nuly a „upper exclusive“ bylo často kritizováno jako nepřirozené, na druhé straně vášnivě obhajováno některými počítačovými teoretiky jako jediné správné. Oba způsoby mohou být přirozenější v různých situacích, a oba jsou již pevně zakořeněny v různých oblastech, takže ani jeden zřejmě v dohledné době nezmizí. Je to podobné jako se slipaři a trenkaři (nebo Linuxem a Windows).

Přiřazení where a forall

Kromě jednoduchého přiřazení polí existují ještě další dvě formy přiřazení: maskované přiřazení where a indexované přiřazení forall.

Maskované přiřazení where má tvar

 where(maska) přiřazení

Maska je logické pole stejných rozměrů jako levá strana, určující elementy, pro které se přiřazení skutečně provede. Příkaz má i blokovou formu, ve které lze použít větve elsewhere, fungující analogicky jako u konstrukce  if.

forall funguje jako přiřazení s explicitními indexy.

forall(index=triplet,...,maska) přiřazení

V přiřazení  musí být levá strana závislá na indexech, pravá může být rovněž. Příkaz se provede tak, jako by se pro všechny aktivní kombinace indexů (tj. ty, pro které je příslušný prvek masky pravdivý) nejprve vyhodnotily levé a pravé strany, a pak se provedla přiřazení. Levé strany se nesmějí překrývat. Jednotlivá vyhodnocení a přiřazení se mohou provést v jakémkoli pořadí.

Příklady:

where(q /= 0) q = 1/q
forall(i=1:size(a,2),j=1:size(a,3)) a(:,i,j) = a(:,j,i)

Je důležité si uvědomit, že forall není „paralelní do“, ačkoliv ho trochu připomíná. Explicitně pracuje s indexy, ale má vlastnosti přiřazení, především povinnost udělat si automaticky potřebná pomocná pole. Nováčci moderního Fortranu často chybně usuzují, že „polní“ aritmetika je má zbavit nutnosti používat cyklus do.

Pole jako argumenty

Je pochopitelně možné předávat pole jako argumenty procedurám. Podle deklarace (a mechanismu předání „za oponou“) rozlišujeme:

• assumed shape. Formální argument deklarujeme s atributem dimension, v němž specifikace každé dimenze má tvar [dolní mez]:. dolní mez lze vynechat (zbude jen dvojtečka) – pak se použije 1. Při volání lze jako skutečný argument uvést pole nebo jeho sekci. Rozměry pole se předají automaticky „za oponou“. Uvnitř procedury lze pole indexovat pomocí nových dolních mezí (to poskytuje abstrakci). Zatímco uvnitř procedury, která pole deklarovala, jsou meze klíčové, při předávání jiným procedurám jsou obvykle relevantní pouze rozměry pole. Proto si procedura specifikuje nové dolní meze a horní meze se dopočítají. Za oponou se tato pole obvykle předají pomocí deskriptoru (struktury obsahující rozměry a rozložení v paměti).
• explicit shape. Formální argument deklarujeme stejně jako automatické pole, tedy uvedeme obě meze ve všech dimenzích. Je tedy potřeba předat rozměry pole do procedury zvlášť (např. jako argumenty). Výhodou je typicky větší rychlost přístupu k poli, protože tato pole jsou udržována souvislá. Skutečný argument nemusí mít stejný rank jako formální, ale musí mít alespoň tolik prvků. Za oponou se předávají adresou prvního prvku. Je-li skutečný argument nesouvislá sekce, udělá se lokální souvislá kopie.
• assumed size. Stejné jako explicit shape, jen v posledním rozměru se místo horní meze použije hvězdička. Tato poslední horní mez je nedefinovaná a nesmí se nikde použít, a to ani implicitně. Lze použít jen ty prvky formálního argumentu, které odpovídají prvkům skutečného argumentu (opět, tyto informace je potřeba do procedury předat jinak). Předávají se stejně jako explicit shape.

Prvotní Fortran umožňoval jen předávání explicit shape, ve Fortranu 77 přibyl způsob assumed size a ve Fortranu 90 assumed shape. Nejpohodlnější je způsob assumed shape, kdy se rozměry předávají automaticky, a lze pracovat i s nesouvislou sekcí bez lokální kopie. Na druhou stranu, protože se zde rozměry a rozložení v paměti určují za běhu, může mít tento styl negativní dopad na rychlost – překladač musí počítat s tím, že pole může být nesouvislé v paměti a může mít libovolnou délku. Naproti tomu explicit shape garantuje rychlý přístup k poli v průběhu procedury, neboť je souvislé a část rozměrů může být známa už při překladu. Navíc se nemusí dodržet rank. Za to ovšem platíme tím, že informace o rozměrech si musíme do procedury předat sami, a pokud je skutečný argument nesouvislý, překladač při volání procedury musí udělat lokální souvislou kopii (zpravidla na stacku). Pro assumed size platí totéž co pro explicit shape.

Vestavěné funkce

Fortran nabízí slušnou kolekci vestavěných funkcí pro práci s poli, např. spread (obdoba repmat v Matlabu), reshape (změna tvaru pole), product (součin prvků), cshift (cyklický posun) a další. Detailní výčet a popis by byl na dlouho, takže zde odkážeme na manuály.

Elementální procedury

Jedním z velikých přínosů Fortranu 95 je možnost definovat si tzv. elementální funkce a podprogramy, které se chovají stejně jako vestavěné operace a matematické funkce. Např.

elemental function erf(x) result(y)
real,intent(in):: x
real:: y
pomocné proměnné, výpočet
y = výsledek
end function

Splňuje-li funkce určité požadavky (nemá vedlejší efekty, jen skalární vstupní argumenty a skalární výsledek), lze ji označit magickým slovem elemental a pak je elementální, tj. jako argumenty lze používat skaláry i pole (viz výše). Různé speciální matematické funkce tyto požadavky obvykle splňují, takže se často deklarují jako elementální, aby se daly pohodlněji používat. Lze definovat i elementální podprogramy – ty mohou mít i výstupní skalární argumenty, ale nelze je samozřejmě použít ve výrazech.

Poznamenejme, že procedury lze také označit jako pure, což znamená bez vedlejších efektů. Ty nelze používat jako elementální, ale pure funkce lze použít v příkazu  forall.

Historické zajímavosti

V dřívějších Fortranech bylo možné používat reálné indexy polí, které se automaticky zaokrouhlovaly. Tato obskurní vlastnost byla zrušena Fortranem 95 (poté, co byla uznána jako prakticky nepoužívaná). Stejný osud potkal i smyčky do s reálnou kontrolní proměnnou.

Hlavním důvodem, proč se ve Fortranu indexuje v kulatých závorkách (stejně jako volání funkce) byl zřejmě především omezený počet znaků dřívějších počítačů (ano, Fortran je starší než ASCII). Také to umožňovalo nahradit funkci s celočíselným argumentem tabulkou předpočítaných hodnot. Z matematického hlediska je to vlastně logické, protože indexování je opravdu jen funkce (zobrazení). Programátorům, z jejichž pohledu jsou indexování pole a volání funkce fundamentálně odlišné věci, to ale často připadá podivné.

Předávání polí stylem assumed size bylo ve Fortranu 95 označeno jako zastaralé. Později však standardizační komise rozhodnutí přehodnotila (mj. i z podnětů uživatelů) a ve Fortranu 2003 už tomu tak není. Assumed size totiž v mnoha případech nejde dobře nahradit zbylými dvěma způsoby.

Příště

Příště si povíme více o typech, mechanismu druhů (kind), odvozených typech (strukturách) a generických rozhraních. Povíme si také o ukazatelích ( pointer). Ukážeme si i práci s řetězci, a zbude-li prostor, podíváme se na některé věci z Fortranu 2003.