Univerzita Jana Evangelisty Purkyně v Ústí nad Labem Přírodovědecká fakulta

(1)

Univerzita Jana Evangelisty Purkyně v Ústí nad Labem Přírodovědecká fakulta

PARALELNÍ PROGRAMOVÁNÍ S APLIKACEMI Martin Lísal

2007

(2)

Studijní opora je určena studentům, kteří jsou zběhlí v programování v jazyce FORTRAN či v jazyce C nebo C++. Nejdříve je proveden úvod do problematiky paralelního programování.

Jsou vysvětleny důvody paralelizace a její efektivita. Dále je uveden přehled úloh a architek- tury počítačů z hlediska paralelizace a přehled prostředků pro paralelizaci numerických výpočtů.

Vlastní část studijní opory se týká úvodu do jedné z nejpoužívanějších technik paralelizace, tzv.

„Message Passing Interfaceÿ (MPI) knihovny. Nejdříve jsou na jednoduchých programech v jazyce FORTRAN 90 demonstrovány základní MPI příkazy. Poté následuje přehled nejpoužíva- nějších MPI příkazů doplněný ilustračními programy. V závěrečné části jsou uvedeny konkrétní příklady paralelizace tří počítačových metod: paralelizace metody paralelního temperingu pro pohyb částice v jednorozměrném silovém poli a při modelování polymerů, paralelizace moleku- lární dynamiky a paralelizace Monte Carlo metody na příkladu Lennardovy-Jonesovy tekutiny.

Recenzovali: Ing. Zdeněk Wágner, CSc.

Ing. Alexandr Babič, Ph.D.

c Martin Lísal, 2007

ISBN 978-80-7044-902-8

(3)

Obsah

1 ÚVOD 5

1.1 Co je to paralelní počítání? . . . . 5

1.2 Proč potřebujeme paralelní počítače? . . . . 5

1.3 Problémy spojené s vývojem paralelního počítání . . . . 6

1.4 Kdy se vyplatí paralelizace? . . . . 6

1.5 Závisí způsob paralelizace výpočtů na architektuře paralelních počítačů? . 7 1.6 Rozdělení paralelních úloh z hlediska jejich spolupráce během výpočtu . . . 10

1.7 SPMD úlohy a strategie paralelizace . . . . 11

2 MPI 13 2.1 Co je to MPI? . . . . 13

2.2 Vytvoření prostředí pro paralelní počítání . . . . 13

2.3 Program typu „Hello World” . . . . 14

2.4 Argumenty příkazů MPI SEND a MPI RECV . . . . 16

2.5 Více o MPI SEND a MPI RECV . . . . 17

2.5.1 Numerická integrace lichoběžníkovou metodou . . . . 17

2.5.2 Paralelní program pro numerickou integraci lichoběžníkovou metodou 18 2.6 Vstup a výstup v paralelních programech . . . . 20

2.6.1 Vstup z terminálu s použitím MPI SEND a MPI RECV . . . . 21

2.6.2 Vstup z terminálu s použitím MPI BCAST . . . . 23

2.6.3 Vstup z terminálu s použitím MPI PACK a MPI UNPACK . . . . 24

2.6.4 Vstup ze souboru . . . . 26

2.6.5 Srovnání jednotlivých metod . . . . 27

2.7 Příkazy MPI REDUCE a MPI ALLREDUCE . . . . 27

2.8 Často používané MPI příkazy . . . . 30

2.8.1 Příkazy pro vytvoření a správu paralelního prostředí . . . . 30

2.8.2 Příkazy pro kolektivní komunikaci . . . . 31

2.8.3 Příkazy pro operace na proměnných distribuovaných na jednotlivých procesech . . . . 39

3 APLIKACE 43 3.1 Paralelní „tempering” . . . . 43

3.1.1 Částice v jednorozměrném silovém poli . . . . 43

3.1.2 Monte Carlo metoda . . . . 44

3.1.3 Metoda paralelního „temperingu” . . . . 49

3.2 Paralelní „tempering” při modelování polymerů . . . . 57

3.2.1 Úvod . . . . 57

3.2.2 Konformace polymeru . . . . 58

3.2.3 Model polymeru . . . . 59

3.2.4 Rosenblutův Monte Carlo algoritmus . . . . 60

3.2.5 Rosenblutovo Monte Carlo a paralelní „tempering” . . . . 61

(4)

3.3 Paralelní molekulární dynamika . . . . 82

3.3.1 Úvod . . . . 82

3.3.2 Paralelizace . . . . 82

3.4 Paralelní Monte Carlo . . . . 93

3.4.1 Úvod . . . . 93

3.4.2 Hybridní Monte Carlo metoda . . . . 94

3.4.3 Paralelizace . . . . 95

(5)

1 ÚVOD

1.1 Co je to paralelní počítání?

Paralelní počítání je počítání na paralelních počítačích či jinak řečeno využití více než jednoho procesoru při výpočtu (viz obr. 1).

(a) Jednoprocesorový poþítaþ

(b) Paralelní poþítaþ typu výpoþtový cluster

CPU memory CPU

memory

CPU memory

Obrázek 1: Schema (a) jednoprocesorového a (b) paralelního počítače typu výpočtový cluster; paměť (memory), procesor (CPU).

1.2 Proč potřebujeme paralelní počítače?

• Od konce 2. světové války dochází k bouřlivému rozvoji počítačových simulací (com-

putational science). Počítačové simulace řeší problémy, které teoretická věda (the-

oretical science) nedokáže vyřešit analyticky či problémy, které jsou obtížné nebo

nebezpečné pro experimentální vědu (experimental science). Jako příklad uveďme

proudění okolo trupu letadla, které teoretická věda dokáže popsat parciálními dife-

renciálními rovnicemi; ty ale nedokáže analyticky řešit. Dalším příkladem je určení

vlastností plazmy. Jelikož plazma existuje při vysokých teplotách, je experimentální

měření většiny jejích vlastností obtížné či přímo nemožné. Počítačové simulace vy-

žadují náročné numerické výpočty, které se v současnosti neobejdou bez nasazení

paralelizace.

(6)

• Vývoj procesorů naráží (či v budoucnu narazí) na svůj materiálový limit t.j. vývoj procesorů nepůjde zvyšovat do nekonečna. (To pull a bigger wagon, it is easier to add more oxen than to grow a gigantic ox.)

• Cena nejvýkonějších (advanced) procesorů na trhu obvykle roste rychleji než jejich výkon. (Large oxen are expensive.)

1.3 Problémy spojené s vývojem paralelního počítání

• Hardware: v současné době lze konstruovat paralelní počítače mající společnou (sdí- lenou) paměť pro maximálně asi 100 procesorů. Trendem je tedy spojení jednoproce- sorových či několikaprocesorových počítačů do clusteru pomocí rychlé komunikační sítě tzv. switche. Pozornost v oblasti hardwaru se proto převážně zaměřuje na vývoj switchů, které umožňují rychlost meziprocesorové komunikace srovnatelnou s rych- lostí komunikace mezi procesorem a vnitřní pamětí počítače.

• Algoritmy: bohužel se ukazuje, že některé algoritmy používané na jednoprocesorových počítačích nelze dobře paralelizovat. To vede k nutnosti vytvářet zcela nové paralelní algoritmy.

• Software: existují prostředky pro automatickou paralelizaci jakými jsou např. para- lelní kompilátory [High Performance Fortran (HPF) Standard] či paralelní numerické knihovny. Platí (a pravděpodobně dlouho bude platit), že nejlepší paralelizace je ta, kterou si uděláme sami.

1.4 Kdy se vyplatí paralelizace?

Zrychlení výpočtu paralelizací lze odhadnout pomocí Amdahlova pravidla:

100 − P + P

n (1)

kde P je část programu, kterou lze paralelizovat a n je počet použitých procesorů. Amdah- lovo pravidlo nám říká, kolik procent výpočtového času bude trvat paralelizovaný program oproti stejnému programu spuštěnému na jednoprocesorovém počítači.

Uveďme si jednoduchý ilustrativní příklad. Představme si, že máme sečíst 1 000 000 000 prvků vektoru a. Načtení prvků trvá 8 % celkového výpočtového času, vlastní výpočet součtu trvá 90 % celkového výpočtového času a tisk trvá 2 % celkového výpočtového času.

Je zřejmé, že jen výpočet součtu lze paralelizovat t.j. P = 90 % (viz. obr. 2). V případě

použití n = 10 procesorů dostaneme z Amdahlova pravidla výsledek 19 % t.j. zrychlení

výpočtu paralelizací přibližně 5krát. Z Amdahlova pravidla a uvedeného příkladu plynou

(7)

dva závěry: (i) neplatí přímá úměra mezi počtem procesorů a urychlením výpočtu („užití 10 procesorů neznamená urychlení výpočtu 10krát”) a (ii) pro limitní případ n → ∞ je urychlení výpočtu konečné a rovno 100/(100 − P ).

(a) Sériový program

S1

P

S2

(b) Paralelní program

a.out

S1

S2

S1

S2

S1

S2

S1

S2

P P P P

Obrázek 2: Schéma (a) sériového a (b) paralelního programu (běžícího na 4 procesorech). S

₁

a S

₂

jsou části programu, které nelze paralelizovat a P je paralelizovatelná část programu.

1.5 Závisí způsob paralelizace výpočtů na architektuře paralel- ních počítačů?

Způsob paralelizace výpočtů bohužel závisí na architektuře paralelních počítačů. V sou- časnosti se nejvíce používají dva typy architektur: (i) paralelní počítače se sdílenou pamětí (shared memory) a (ii) paralelní počítače s distribuovanou pamětí (distributed memory);

viz obr. 3.

(8)

(a) Paralelní poþítaþ se sdílenou pamČtí (shared memory)

(b) Paralelní poþítaþ s distribuovanou pamČtí

(distributed memory)

(shared) memory

CPU CPU CPU CPU memory

CPU

memory

CPU

memory

CPU

memory

CPU switch ~ komunikace

Obrázek 3: Schema paralelního počítače (a) se sdílenou pamětí (shared memory) a (b) s distribuovanou pamětí (distributed memory).

Z hlediska způsobu paralelizace je výhodnější architektura se sdílenou pamětí, neboť sdílení a výměna informací mezi procesory během výpočtu se děje přes sdílenou (společ- nou) paměť. Paralelizace na „shared-memory” počítačích se provádí pomocí tzv. OpenMP.

OpenMP je seznam direktiv, které se vkládají na místa, které chceme paralelizovat. Jinými slovy říkáme počítači „kde a co” paralelizovat. Uveďme si jako příklad paralelizaci cyklu (do) pomocí OpenMP:

!$OMP PARALLEL DO do i = 2, n

b(i) = (a(i) + a(i-1)) / 2.0 end do

!$OMP END PARALLEL DO

Direktiva před začátek cyklu, !$OMP PARALLEL DO, říká: „paralelizuj cyklus”, zatímco direktiva na konci cyklu, !$OMP END PARALLEL DO, říká: „ukonči paralelizaci cyklu”.

¹

Paralelní počítače se sdílenou pamětí jsou oproti paralelním počítačům s distribuovanou pamětí dražší a počet procesorů sdílející jednu paměť je v současnosti omezen maximálně asi na 100. Proto se běžněji setkáváme s „distributed-memory” počítači. „Distributed- memory” počítače mají fyzicky oddělené paměti a komunikace (sdílení dat) mezi procesory

1

Jelikož direktivy OpenMP začínají ”!” a ten se ve FORTRANu90 interpretuje jako začátek komentářo-

vého řádku, lze programy s direktivami OpenMP spouštět bez problému na jednoprocesorových počítačích.

(9)

zajišťuje switch. Paralelizace na „distributed-memory” počítačích je obtížnější než para- lelizace na „shared-memory” počítačích, neboť musíme zajistit a pracovat s předáváním informací (message passing) mezi procesory během výpočtu. Na paralelních počítačích s distribuovanou pamětí se používají dva systémy: (i) Parallel Virtual Machine (PVM) vhodný pro heterogenní clustery a (ii) Message Passing Interface (MPI) vhodný pro homo- genní clustery. PVM a MPI je knihovna příkazů pro FORTRAN, C a C++, jež umožňují předávání informací mezi procesory během výpočtu.

Komunikace představuje časovou ztrátu, o kterou se sníží zrychlení výpočtu paralelizací dané Amdahlovým pravidlem. Čas komunikace (overhead) jako funkce množství předávané informace je znázorněn na obr. 4. Z obr. 4 je patrno, že čas komunikace závisí na latenci (času, který potřebuje procesor a switch k „navázání komunikace”) a dále, že „overhead”

roste s množstvím předávané informace mezi procesory. Z obr. 4 plynou následující zásady, které bychom měli dodržovat při paralelizaci na počítačích s distribuovanou pamětí: (i) přenášet co nejméně informací a (ii) komunikovat co nejméně.

množství pĜenášené informace

þas komunikace (overhead)

latence

Obrázek 4: Čas komunikace (overhead) versus množství předávané informace mezi proce-

sory.

(10)

1.6 Rozdělení paralelních úloh z hlediska jejich spolupráce bě- hem výpočtu

Podle spolupráce během výpočtu můžeme rozdělit paralelní úlohy na

• MPMD (Multiple Program Multiple Data) úlohy

• SPMD (Single Program Multiple Data) úlohy

MPMD úlohy lze dále rozdělit na typ

• „Master/Worker” (obr. 5a)

• „Coupled Multiple Analysis” (obr. 5b)

Úlohy typu „Master/Worker” jsou charakteristické existencí řídícího procesu (master, a.out v obr. 5a), který řídí ostatní procesy (worker, b.out v obr. 5a) a nechává si od nich počítat mezivýsledky. Úlohy typu „Coupled Multiple Analysis” jsou charakteristické použitím nezávislých procesů (a.out, b.out a c.out v obr. 5b), které řeší různé aspekty určitého problému. Např. při návrhu trupu letadla by program a.out řešil proudění, b.out by řešil pevnostní analýzu a c.out analýzu teplotní.

SPMD úlohy (obr. 5c) používají stejnou kopii programu (single program) na všech

procesorech. Procesy ale během výpočtu zpracovávají různá data (multiple data). SPMD

úlohy jsou typické pro numerické výpočty a my se SPMD úlohami budeme výhradně za-

bývat v rámci tohoto předmětu.

(11)

(a) MPMD úloha typu "Master/Worker"

(b) MPMD úloha typu

"Coupled Multiple Analysis"

a.out b.out c.out

b.out a.out

(c) SPMD úloha a.out

Obrázek 5: Paralelní úlohy z hlediska jejich spolupráce během výpočtu. (a) MPMD (Multiple Program Multiple Data) úloha typu „Master/Worker”, (b) MPMD úloha typu

„Coupled Multiple Analysis” a (c) SPMD (Single Program Multiple Data) úloha.

1.7 SPMD úlohy a strategie paralelizace

Uveďme si jednoduchý příklad paralelizace programu pro součet prvků vektoru, který ozřejmí pojem SPMD úlohy a naznačí strategii paralelizace SPMD úloh. Představme si, že máme vektor a mající 9 000 000 000 prvků. Úkol je sestavit paralelní program pro sou- čet prvků vektoru a ( ^P

9 000 000 000

i=1

a

_i

), přičemž máme k dispozici 3 procesory. Rozdělme si úlohu pomyslně na 3 části: (i) načtení vektoru a, S

₁

, (ii) vlastní výpočet součtu, P , a (iii) tisk součtu, S

₂

. Je zřejmé, že jen vlastní výpočet součtu P lze paralelizovat.

²

Paralelizaci provedeme tak, že každý procesor spočte 1/3 součtu a tu pak pošle na jeden procesor, který provede celkový součet, jež následně vytiskne. Paralelní program pro součet prvků vektoru zapsaný v jazyce FORTRAN může vypadat následovně (v místech programu, jež by měla obsahovat MPI paralelní příkazy, jsou komentáře označeny jako „. . .”):

program soucet implicit none

integer :: i,my_rank,p,n_loc,i_begin,i_end real :: a(9000000000), &

sum_loc,sum_tot

!

2

Načtení dat se nedoporučuje paralelizovat, neboť způsob načítání dat se liší systém od systému.

(12)

open(10,file=’input.dat’) do i=1,9000000000

read(10,*) a(i) end do

close(10)

!

"MPI: zjisti počet procesů ’p’, které se podílí na výpočtu"

"MPI: zjisti moje pořadí ’my_rank’ ve výpočtu; my_rank=0,1,...,p-1"

!

n_loc=9000000000/p ! Předpokládáme dělitelnost beze zbytku :-) i_begin=my_rank*n_loc+1

i_end=(my_rank+1)*n_loc sum_loc=0.0

do i=i_begin,i_end sum_loc=sum_loc+a(i) end do

!

"MPI: zajisti sečtení jednotlivých ’sum_loc’ do ’sum_tot’

a pošli ’sum_tot’ např. na proces mající ’my_rank=0’"

!

if(my_rank == 0) then write(*,*) sum_tot end if

stop "soucet: Konec vypoctu!"

!

end program soucet

(13)

2 MPI

2.1 Co je to MPI?

MPI (Message Passing Interface) je knihovna funkcí a podprogramů, která umožňuje:

• vytvořit prostředí pro paralelní počítání;

• komunikaci mezi procesy během výpočtu např. posílání informací z procesu na proces;

• kolektivní operace např. sečtení mezivýsledků z procesů a uložení celkové součtu na určitý proces;

• topologizaci procesů t.j. např. seskupení určitého počtu procesů za účelem jejich spolupráce v rámci výpočtů.

MPI lze použít ve spojení s programovacími jazyky FORTRAN, C a C++.

2.2 Vytvoření prostředí pro paralelní počítání

Vytvořením prostředí pro paralelní počítání rozumíme:

• inicializaci a ukončení používání MPI;

• definování celkového počtu procesů nprocs, které se účastní výpočtu;

• definování pořadí procesu myrank.

Demonstrujme si vytvoření paralelního prostředí na jednoduchém programu, který vy- píše na obrazovku pořadí jednotlivých procesů.

³

program par_env implicit none

include ’mpif.h’ ! Budeme pouzivat MPI knihovnu

!

integer :: ierr,nprocs,myrank

!

call MPI_INIT(ierr) ! INITializuje MPI

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr)

! Default COMMunikator a celkovy pocet (SIZE) procesu ’nprocs’

!

call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ierr)

! Poradi (RANK) ’myrank’; myrank = 0, 1, ..., nprocs-1

!

print*," Poradi procesu =",myrank

!

call MPI_FINALIZE(ierr) ! Ukonci (FINALIZE) MPI end program par_env

3

Pokud nebude řečeno jinak, pracujeme s SPMD modelem paralelního programu: na všech uzlech se

rozběhne stejná kopie programu (SP, single program), každá kopie programu pracuje s různými daty (MD,

multiple data).

(14)

Komentář k programu:

• Include ’mpif.h’ říká, že budeme používat MPI knihovnu.

• MPI INIT(ierr) inicializuje MPI; ierr je „chybová” proměnná typu INTEGER (ierr = 0 když je vše v pořádku).

• MPI COMM SIZE(MPI COMM WORLD,nprocs,ierr) definuje celkový po- čet procesů nprocs, které se účastní výpočtu.

• MPI COMM RANK(MPI COMM WORLD,myrank,ierr) definuje pořadí procesu myrank.

• MPI FINALIZE(ierr) ukončí MPI.

• V celém programu používáme předdefinovaný (default) tzv. „komunikátor”

MPI COMM WORLD – označení pro skupinu spolupracujících procesů.

2.3 Program typu „Hello World”

Při učení nového programovacího jazyka se pravděpodobně každý setkal s programem, který vypíše na obrazovku „pozdrav” - program typu „Hello World”. Napišme si takovou verzi programu pomocí MPI. MPI „Hello World” program pošle z procesů majících pořadí myrank = 1, 2, . . . , nprocs −1 informaci o svém pořadí na proces s pořadím 0 a proces s pořadím 0 vytiskne tuto informaci na obrazovku.

program hello_world implicit none

include ’mpif.h’ ! Budeme pouzivat MPI knihovnu

!

integer :: STATUS(MPI_STATUS_SIZE)

! Pomocna promenna pouzita v MPI_RECV integer :: ierr,nprocs,myrank,tag,isource,act_proc

!

call MPI_INIT(ierr) ! INITializujeme MPI

!

! Default COMMunikator a celkový pocet (SIZE) procesu ’nprocs’

!

! Poradi (RANK) ’myrank’; myrank = 0, 1, ..., nprocs-1

! tag=1

if(myrank /= 0) then ! Ne_Mastr procesy act_proc=myrank

call MPI_SEND(act_proc,1,MPI_INTEGER,0,tag,MPI_COMM_WORLD,ierr)

! Ne_Mastr procesy SEND info ’act_proc’

else ! Mastr proces do isource=1,nprocs-1

call MPI_RECV(act_proc,1,MPI_INTEGER,isource,tag,MPI_COMM_WORLD, &

STATUS,ierr)

! Mastr proces RECV info ’act_proc’

(15)

print*," Hello World =",act_proc enddo

endif

!

call MPI_FINALIZE(ierr) ! Ukonci (FINALIZE) MPI end program hello_world

Komentář k programu:

• MPI SEND(. . . ) posílá informace.

• MPI RECV(. . . ) zajistí obdržení informace.

• 1. argument v MPI SEND/MPI RECV je posílaná/obdržená informace.

• 2. argument v MPI SEND/MPI RECV je velikost posílané/obdržené informace (1 pro skalár, počet prvků vektoru či matice, atd.).

• 3. argument v MPI SEND/MPI RECV značí MPI typ, který přibližně odpovídá ty- pům v jazycích FORTRAN, C a C++.

• 4. argument v MPI SEND značí, na jaký proces se má zpráva poslat (v našem případě na proces 0).

• 4. argument v MPI RECV značí, z jakých procesů zpráva přichází (v našem případě z procesů 1, 2, . . . , nprocs −1).

• 5. argument v MPI SEND/MPI RECV je proměnná typu INTEGER (0 – 32767), která umožňuje dodatečné (pomocné) označení zprávy.

• Proměnná STATUS v MPI RECV zpětně „informuje” MPI SEND o úspěšném do- ručení zprávy.

⁴

• MPI RECV je až na STATUS zrcadlovým obrazem MPI SEND.

• MPI SEND a MPI RECV jsou vzájemně synchronizované příkazy t.j. běh programu pokračuje, až když se oba příkazy dokončí.

4

Ostatní argumenty v MPI SEND/MPI RECV byly již vysvětleny.

(16)

2.4 Argumenty příkazů MPI SEND a MPI RECV

V programu typu „Hello World” jsme použili příkaz MPI SEND pro posílání informace z uzlu na uzel a příkaz MPI RECV pro přijímání informací z uzlů. Popišme si detailně argumenty těchto příkazů. Tyto dva příkazy lze obecně zapsat jako:

call MPI_SEND(send_message,count,MPI_data_type,dest,tag,comm, &

ierr)

call MPI_RECV(recv_message,count,MPI_data_type,source,tag,comm, &

status,ierr)

kde

• send message je posílaná informace, jejíž typ specifikuje MPI data type a jež může být skalár, vektor či matice.

• recv message je přijímaná informace, jejíž typ specifikuje MPI data type a jež může být skalár, vektor či matice.

• count je typu INTEGER a vyjadřuje velikost posílané resp. přijímané informace;

např. 1 pro skalár či 10 pro matici (2,5).

• MPI data type :

MPI data type FORTRAN

MPI INTEGER INTEGER

MPI REAL REAL

MPI DOUBLE PRECISION DOUBLE PRECISION

MPI COMPLEX COMPLEX

MPI LOGICAL LOGICAL

MPI CHARACTER CHARACTER

MPI PACKED MPI BYTE

• dest je typu INTEGER a vyjadřuje pořadí procesu, na který se send message posílá.

• source je typu INTEGER a vyjadřuje pořadí procesu, z kterého se recv message přijímá.

• tag je proměnná typu INTEGER, která může nabývat hodnot od 0 do 32767 a která umožňuje dodatečné (pomocné) označení posílané resp. přijímané informace.

• comm je označení pro skupinu spolupracujících procesů („komunikátor”);

MPI COMM WORLD je předdefinovaný (default) „komunikátor”.

(17)

• status je vektor typu INTEGER deklarovaný jako integer :: status(MPI_STATUS_SIZE)

status zpětně „informuje” MPI SEND o úspěšném doručení zprávy.

• ierr je „chybová” proměnná typu INTEGER; ierr = 0 když je vše v pořádku.

2.5 Více o MPI SEND a MPI RECV

Ukažme si další použití MPI SEND a MPI RECV na paralelní verzi programu pro numerický výpočet určitého integrálu lichoběžníkovou metodou.

2.5.1 Numerická integrace lichoběžníkovou metodou

Úkolem je spočíst integrál funkce f (x):

I =

Z

b a

f (x) dx (2)

Krok 1: rozdělme interval ha, bi na n stejných dílů:

x

_i

= a + ih i = 0, 1, 2, . . . , n − 1, n (3)

x

₀

= a (4)

x

_n

= b (5)

kde

h = b − a

n (6)

Krok 2: aproximujme funkci f (x) mezi body x

_i

a x

_i+1

přímkou.

Krok 3: body x

_i

, x

_i+1

, f (x

_i

) a f (x

_i+1

) tvoří lichoběžník, jehož obsah je A

_i

= h f (x

i

) + f (x

i+1

)

2 (7)

Krok 4: sečtením obsahů všech lichoběžníků, dostaneme přibližnou hodnotu integrálu (2):

I ≈ h

"

f (a) + f (b)

2 +

n−1

X

i=1

f (x

_i

)

#

(8)

≈

"

f (a)

2 + f (a + h) + f (a + 2h) + . . . + f (b − h) + f (b) 2

#

(9)

(18)

Sériový program pro numerickou integraci lichoběžníkovou metodou může vypadat ná- sledovně:

program ser_lich_int implicit none

!

integer :: i,n

real :: f,a,b,h,x,int_sum,Int

!

print *,"Dolni mez intervalu a:"

read *,a

print *,"Horni mez intervalu b > a:"

read *,b

print *,"Zadejte deleni intervalu n > 0:"

read *,n

!

h=(b-a)/REAL(n) ! "krok"

int_sum=(f(a)+f(b))/2.0 ! [f(a)+f(b)]/2 do i=1,n-1

x=a+REAL(i)*h ! x=a+i*h

int_sum=int_sum+f(x) ! pomocna suma end do

Int=h*int_sum ! vysledny integral

!

print *,"Integral =",Int end program ser_lich_int

!

function f(x) implicit none

!

real :: f,x f=x*x end function f

2.5.2 Paralelní program pro numerickou integraci lichoběžníkovou metodou Paralelní verze programu pro numerickou integraci lichoběžníkovou metodou

1. rozdělí interval ha, bi na podintervaly ha

_i

, b

_i

i dle celkového počtu procesů nprocs, 2. pro jednotlivé podintervaly spočtou jednotlivé procesy integrály

Z

bi

ai

f (x) dx pomocí lichoběžníkové metody

3. a tyto integrály se pošlou na proces mající pořadí 0, kde se sečtou a proces mající pořadí 0 součet (celkový integrál) vytiskne.

Paralelní program pro numerickou integraci lichoběžníkovou metodou může vypadat následovně:

^{5 6}

5

Abychom se vyhnuli problematice načítání vstupních dat v paralelních programech, zadejme vstupní hodnoty a, b a n přímo do programu.

6

Dále předpokládáme, že n je beze zbytku dělitelné nprocs.

(19)

!

! Trapezoidal Rule, MPI version 1.

!

! a, b, and n are in DATA statement

!

! Algorithm:

! 1. Each process calculates "its" interval of integration.

! 2. Each process estimates the integral of f(x) over its interval

! using the trapezoidal rule.

! 3a. Each process /= 0 sends its integral to 0.

! 3b. Process 0 sums the calculations received from the individual

! processes and prints the result.

!

program par_lich_int1 implicit none

include ’mpif.h’ ! preprocessor directive

!

integer :: nprocs, & ! # of processes myrank, & ! my process rank

source, & ! process sending integer dest, & ! process receiving integer tag,ierr

integer :: status(MPI_STATUS_SIZE) integer :: n, & ! # of trapezoids

local_n ! # of trapezoids for a processor

!

real :: a, & ! left endpoint b, & ! right endpoint h, & ! trapezoid base length

local_a, & ! left endpoint for a processor local_b, & ! right endpoint for a processor integral, & ! integral over a processor total ! total integral

!

data a,b,n /0.0,1.0,1024/ ! data for integration, n must be even data dest,tag /0,50/

! ! start up MPI ! call MPI_INIT(ierr)

! ! find out how many processes are being used ! call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr)

! ! get my process rank !

!

if(mod(n,nprocs) /= 0) stop "par_lich_int1: Wrong n/nprocs ratio!"

h=(b-a)/REAL(n) ! h is the same for all processes local_n=n/nprocs ! # of trapezoids

! ! length of interval of integration for a process ! local_a=a+REAL(myrank)*REAL(local_n)*h

local_b=local_a+REAL(local_n)*h

! ! calculate integral on a process !

call Trap(local_a,local_b,local_n,h,integral)

! ! add up the integrals calculated by each process ! if(myrank /= 0) then

call MPI_SEND(integral,1,MPI_REAL,dest,tag, &

MPI_COMM_WORLD,ierr) else

total=integral do source=1,nprocs-1

call MPI_RECV(integral,1,MPI_REAL,source,tag, &

MPI_COMM_WORLD,status,ierr) total=total+integral

enddo endif

! ! print the results ! if(myrank == 0) then

(20)

write(*,’(1x,a,e13.5)’) " Integral =",total endif

! ! shut down MPI ! call MPI_FINALIZE(ierr) end program par_lich_int1

! ! Subroutine trapezoid for a processor !

subroutine Trap(local_a,local_b,local_n,h,integral) implicit none

!

integer :: local_n,i

real :: local_a,local_b,h,integral,x,f

!

integral=(f(local_a)+f(local_b))/2.0 x=local_a

do i=1,local_n-1 x=x+h

integral=integral+f(x) enddo

integral=integral*h end subroutine Trap

! ! Function for integration ! function f(x)

implicit none

!

2.6 Vstup a výstup v paralelních programech

MPI standard neobsahuje speciální příkazy pro načítání dat či jejich výstup na terminál nebo do souboru. Vstup a výstup v paralelních programech se proto musí řešit pomocí příkazů programovacího jazyka read, write či print a MPI příkazů pro komunikaci mezi uzly.

Je nutno si uvědomit, že např. způsob načtení dat použitý v sériovém programu pro numerickou integraci lichoběžníkovou metodou:

print *,"Dolni mez intervalu a:"

read *,a

print *,"Horni mez intervalu b > a:"

read *,b

print *,"Zadejte deleni intervalu n > 0:"

read *,n

může v paralelním programu fungovat různě v závislosti na systému a architektuře počítače.

Při podobném použití příkazu read v paralelních programech může dojít k následujícím situacím:

• Pokud jen jeden proces „umí” číst z terminálu, data se načtou jen na tento proces a

ne na ostatní procesy.

(21)

• Pokud „umí” všechny procesy číst, pak může nastat případ zadávat data tolikrát, kolik procesorů používáme.

⁷

Situace se dále komplikuje, pokud načítáme data ze souboru a např. více procesů musí číst současně data z jednoho souboru. Obdobné problémy a nejednoznačnosti jsou spojeny s výstupem na terminál či do souborů.

Systémy a architektury počítačů vždy minimálně umožňují následující:

• Alespoň jeden proces dokáže číst z terminálu a alespoň jeden proces dokáže zapisovat na terminál.

• Každý proces dokáže číst ze souboru či do souboru zapisovat, pokud ostatní procesy z tohoto souboru nečtou či do něj nezapisují.

• Procesy dokáží číst současně z různých souborů či do různých souborů zapisovat.

Na základě těchto „minimálních” schopností systémů a architektur počítačů se doporu- čuje používat následující dvě zásady při řešení vstupu a výstupu do paralelních programů přes terminál:

1. Načíst data jen jedním procesem (např. procesem s myrank = 0) a následně rozeslat data z tohoto procesu na ostatní procesy.

2. Poslání výsledků výpočtů z jednotlivých procesů na jeden určitý proces (např. proces s myrank = 0) a tisk výsledků z tohoto procesu na terminál.

V následujícím si ukažme různá řešení vstupů dat do paralelních programů z terminálu a souboru na příkladu paralelního programu pro numerickou integraci lichoběžníkovou metodou. Abychom se vyhnuli opakovanému psaní paralelního programu pro numerickou integraci lichoběžníkovou metodou, nahradíme příkaz

data a,b,n /0.0,1.0,1024/ ! data for integration, n must be even

podprogramem, který provede načtení dat. Zbytek programu se nezmění.

2.6.1 Vstup z terminálu s použitím MPI SEND a MPI RECV

V první verzi vstupu z terminálu načte proces s pořadím myrank = 0 data a tato data pošle s použitím příkazů MPI SEND a MPI RECV na ostatní procesy. Příslušný podpro- gram může vypadat následovně:

!

! Algorithm:

! 1. Process 0 reads data.

7

Dále není zaručeno, že čtení proběhne podle pořadí procesů.

(22)

! 2. Process 0 sends data to other processes.

!

subroutine Get_Data1(a,b,n,myrank,nprocs) implicit none

!

integer :: myrank, & ! my process rank nprocs, & ! # of processes

integer :: n ! # of trapezoids integer :: status(MPI_STATUS_SIZE)

!

real :: a, & ! left endpoint

b ! right endpoint

!

if(myrank == 0) then print *," a:"

read*,a print *," b:"

read*,b print *," n:"

read*,n

!

do dest=1,nprocs-1 tag=0

call MPI_SEND(a,1,MPI_REAL,dest,tag,MPI_COMM_WORLD,ierr) tag=1

call MPI_SEND(b,1,MPI_REAL,dest,tag,MPI_COMM_WORLD,ierr) tag=2

call MPI_SEND(n,1,MPI_INTEGER,dest,tag,MPI_COMM_WORLD,ierr) enddo

else tag=0

call MPI_RECV(a,1,MPI_REAL,0,tag,MPI_COMM_WORLD, &

status,ierr) tag=1

call MPI_RECV(b,1,MPI_REAL,0,tag,MPI_COMM_WORLD, &

status,ierr) tag=2

call MPI_RECV(n,1,MPI_INTEGER,0,tag,MPI_COMM_WORLD, &

status,ierr) endif

end subroutine Get_Data1

Tato verze programu volá několikrát příkazy MPI SEND a MPI RECV. Někdy proto bývá výhodnější načtená data nejdříve uložit do pomocného vektoru a ten poslat s použitím příkazů MPI SEND a MPI RECV na ostatní procesy. V takovém případě se sníží latentní čas komunikace. Příslušný podprogram může vypadat následovně:

!

! Algorithm:

! 1. Process 0 reads data and save them to a vector.

! 2. Process 0 sends data as vector to other processes.

!

integer :: myrank, & ! my process rank

(23)

nprocs, & ! # of processes

integer :: n ! # of trapezoids integer :: status(MPI_STATUS_SIZE)

!

real :: a, & ! left endpoint b, & ! right endpoint temp_vec(3) ! temporary vector

!

read*,n

!

temp_vec=(/a,b,REAL(n)/)

!

do dest=1,nprocs-1 tag=0

call MPI_SEND(temp_vec,3,MPI_REAL,dest,tag,MPI_COMM_WORLD, &

ierr) enddo

else tag=0

call MPI_RECV(temp_vec,3,MPI_REAL,0,tag,MPI_COMM_WORLD,status, &

ierr) a=temp_vec(1) b=temp_vec(2) n=INT(temp_vec(3)) endif

2.6.2 Vstup z terminálu s použitím MPI BCAST

V druhé verzi vstupu z terminálu načte opět proces s pořadím myrank = 0 data a tato data se rozešlou s použitím příkazu MPI BCAST

⁸

na ostatní procesy. Příslušný podpro- gram může vypadat následovně:

!

! Algorithm:

! 2. Process 0 broadcasts data to other processes.

!

subroutine Get_Data2(a,b,n,myrank) implicit none

!

integer :: myrank, & ! my process rank ierr

integer :: n ! # of trapezoids

!

b ! right endpoint

8

Příkaz MPI BCAST si můžeme představit jako „sloučení” příkazů MPI SEND a MPI RECV do jed-

noho.

(24)

!

read*,n endif

call MPI_BCAST(a,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(b,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) end subroutine Get_Data2

Příkaz MPI BCAST lze obecně zapsat jako:

call MPI_BCAST(broadcast_message,count,MPI_data_type,root,comm,ierr)

kde význam všech argumentů kromě root je zřejmý z předchozích probíraných MPI příkazů.

root je typu INTEGER a vyjadřuje pořadí procesu, z kterého se broadcast message rozesílá.

2.6.3 Vstup z terminálu s použitím MPI PACK a MPI UNPACK

V třetí verzi vstupu z terminálu načte opět proces s pořadím myrank = 0 data. Tato data před rozesláním „zabalí” do pomocné proměnné buffer s použitím příkazu MPI PACK.

Proměnná buffer se pak rozešle na ostatní procesy s pomocí příkazu MPI BCAST, kde se zpětně „rozbalí” s použitím příkazu MPI UNPACK.

⁹

Příslušný podprogram může vypadat následovně:

!

! Use of Pack/Unpack

!

integer :: myrank, & ! my process rank position, & ! in the buffer ierr

!

b ! right endpoint

!

character(len=100) :: buffer

!

9

Užití MPI PACK a MPI UNPACK je výhodné při posílání velkého množství dat různých typů.

(25)

read*,n

!

! Pack the data into buffer.

! Beginning of buffer ’position=0’

!

position=0 ! Position is in/out

call MPI_PACK(a,1,MPI_REAL,buffer,100,position, &

MPI_COMM_WORLD,ierr)

! Position has been incremented;

! the first free location in buffer call MPI_PACK(b,1,MPI_REAL,buffer,100,position, &

! Position has been incremented again call MPI_PACK(n,1,MPI_INTEGER,buffer,100,position, &

! Position has been incremented again

!

! Broadcast contents of buffer

!

call MPI_BCAST(buffer,100,MPI_PACKED,0,MPI_COMM_WORLD, &

ierr) else

call MPI_BCAST(buffer,100,MPI_PACKED,0,MPI_COMM_WORLD, &

ierr)

!

! Unpack the data from buffer

!

position=0

call MPI_UNPACK(buffer,100,position,a,1,MPI_REAL, &

call MPI_UNPACK(buffer,100,position,b,1,MPI_REAL, &

call MPI_UNPACK(buffer,100,position,n,1,MPI_INTEGER, &

MPI_COMM_WORLD,ierr) endif

Příkazy MPI PACK a MPI UNPACK lze obecně zapsat jako:

call MPI_PACK(pack_data,count,MPI_data_type,buffer,size,position, &

comm,ierr)

call MPI_UNPACK(buffer,size,position,unpack_data,count,MPI_data_type, &

comm,ierr)

kde význam všech argumentů kromě buffer, size a position je zřejmý z předchozích pro- bíraných MPI příkazů. Všiměte si, že proměnné „zabalené” pomocí příkazu MPI PACK (buffer v našem příkladu) jsou v MPI příkazech typu MPI PACKED.

• buffer je pomocná proměnná, jež musí být deklarovaná jako typ CHARACTER.

(26)

• size je typu INTEGER a vyjadřuje velikost („délku”) pomocné proměnné buffer.

¹⁰

• position je typu INTEGER a vyjadřuje pozici v buffer, od které se pack data ukládají či unpack data jsou uloženy. Před prvním použitím MPI PACK či MPI UNPACK je nutno position vynulovat a pak je position automaticky aktualizována voláním MPI PACK či MPI UNPACK.

2.6.4 Vstup ze souboru

Předpokládejme, že vstupní data a, b a n jsou v souboru input.dat. Aby došlo ke správ- nému načtení dat jednotlivými procesy, je nutno zajistit načítání dat z tohoto souboru postupně. To lze docílit použitím příkazu MPI BARRIER. Příslušný podprogram může vypadat následovně:

!

! Reading from file

!

integer :: myrank, & ! my process rank nprocs, & ! # of processes irank,ierr

!

b ! right endpoint

!

do irank=0,nprocs-1 if(irank == myrank) then

open(10,file="input.dat",status="old") read(10,*) a

read(10,*) b read(10,*) n close(10) endif

call MPI_BARRIER(MPI_COMM_WORLD,ierr) enddo

Příkaz MPI BARRIER funguje tak, že program pokračuje až poté, když všechny procesy zavolají tento příkaz.

Obecný tvar příkazu MPI BARRIER je:

call MPI_BARRIER(comm,ierr)

kde význam argumentů je zřejmý z předchozích probíraných MPI příkazů.

10

size musí být větší než součet velikostí („délek”) „balených” dat.

(27)

2.6.5 Srovnání jednotlivých metod

Použití výše zmíněných metod závisí na množství a typech rozesílaných dat:

• V případě rozesílání malého množství dat různých typů, několikanásobné použití MPI SEND a MPI RECV či MPI BCAST je nejjednodušším řešením.

• V případě rozesílání většího množství dat podobných typů (např. REAL, INTE- GER), uložení dat do pomocného vektoru a použití MPI SEND a MPI RECV či MPI BCAST se jeví jako nejvýhodnější řešení.

• V případě rozesílání velkého množství dat podobných typů (např. REAL, INTE- GER) či většího množství dat různých typů (např. REAL, INTEGER, COMPLEX, CHARACTER) je užití MPI PACK a MPI UNPACK nejlepším řešením.

2.7 Příkazy MPI REDUCE a MPI ALLREDUCE

Příkazy MPI REDUCE a MPI ALLREDUCE umožňují provést určité typy operací (např. sčítání či násobení) na proměnných („mezivýsledcích”), které se nachází na jednot- livých procesech. Zároveň umožňují výsledek těchto operací buď poslat na určitý proces (MPI REDUCE) či rozeslat na všechny procesy (MPI ALLREDUCE). Ukažme si opět použití těchto příkazů v paralelním programu pro numerickou integraci lichoběžníkovou metodou. Připomeňme si, že paralelní program pro numerickou integraci lichoběžníkovou metodou: (i) rozdělil celkový interval na podintervaly podle celkového počtu procesů; (ii) jednotlivé procesy spočetly integrály pro jednotlivé podintervaly pomocí lichoběžníkové metody a (iii) tyto integrály se poslaly pomocí MPI SEND a MPI RECV na proces ma- jící pořadí 0, kde se sečetly a výsledek se vytiskl. Bod (iii) nyní nahradíme příkazem MPI REDUCE a modifikovaná verze paralelního programu pro numerickou integraci li- choběžníkovou metodou může vypadat následovně:

!

! Trapezoidal Rule, MPI version 6.

!

! Algorithm:

! 1. Each process calculates "its" interval of integration.

! 2. Each process estimates the integral of f(x) over its interval

! using the trapezoidal rule.

! 3a. MPI_REDUCE sums the integrals on 0.

! 3b. Process 0 prints the result.

!

program par_lich_int6 implicit none

!

integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr

integer :: n, & ! # of trapezoids

local_n ! # of trapezoids for a processor

!

(28)

real :: a, & ! left endpoint b, & ! right endpoint h, & ! trapezoid base length

local_a, & ! left endpoint for a processor local_b, & ! right endpoint for a processor integral, & ! integral over a processor total ! total integral

!

! start up MPI

!

call MPI_INIT(ierr)

!

! find out how many processes are being used

!

! get my process rank

!

! get data

!

call Get_Data2(a,b,n,myrank)

!

if(mod(n,nprocs) /= 0) stop "par_lich_int6: Wrong n/nprocs ratio!"

h=(b-a)/REAL(n) ! h is the same for all processes local_n=n/nprocs ! # of trapezoids

!

! length of interval of integration for a process

!

local_a=a+REAL(myrank)*REAL(local_n)*h local_b=local_a+REAL(local_n)*h

!

! calculate integral on a process

!

call Trap(local_a,local_b,local_n,h,integral)

!

! add up the integrals calculated by each process using MPI_REDUCE

!

call MPI_REDUCE(integral,total,1,MPI_REAL,MPI_SUM,0, &

!

! print the results

!

if(myrank == 0) then

write(*,’(1x,a,e13.5)’) " Integral =",total endif

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program par_lich_int6

!

! Subroutine trapezoid for a processor

!

subroutine Trap(local_a,local_b,local_n,h,integral) implicit none

!

integer :: local_n,i

real :: local_a,local_b,h,integral,x,f

!

integral=(f(local_a)+f(local_b))/2.0 x=local_a

do i=1,local_n-1 x=x+h

(29)

integral=integral+f(x) enddo

integral=integral*h end subroutine Trap

!

! Algorithm:

! 2. Process 0 broadcasts data to other processes.

!

integer :: myrank, & ! my process rank ierr

!

b ! right endpoint

!

read*,n endif

call MPI_BCAST(a,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(b,1,MPI_REAL,0,MPI_COMM_WORLD,ierr) call MPI_BCAST(n,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) end subroutine Get_Data2

!

! Function for integration

!

function f(x) implicit none

!

Argument MPI SUM v příkazu MPI REDUCE říká, že se má provést součet proměnných integral z jednotlivých procesů do proměnné total. Šestý argument 0 v příkazu MPI REDUCE říká, že se výsledek total (součet proměnných integral ) má poslat na proces s pořadím 0.

Obecný tvar příkazu MPI REDUCE lze obecně zapsat jako:

call MPI_REDUCE(operand,result,count,MPI_data_type,operation,root, &

comm,ierr)

kde význam argumentů kromě operand, result, operation a root je zřejmý z předchozích probíraných MPI příkazů.

• operand je proměnná, s kterou MPI REDUCE provádí určité typy operací.

(30)

• result je proměnná typu MPI data type, která obsahuje výsledek určitého typu ope- race na proměnných operand.

• operation říká, jaký typ operace se provádí:

¹¹

operation Význam

MPI SUM součet, ^P

MPI PROD násobení, ^Q

MPI MAX maximum

MPI MIN minimum

MPI MAXLOC maximum a pořadí procesu, kde se maximum nachází MPI MINLOC minimum a pořadí procesu, kde se minimum nachází

• root je typu INTEGER a vyjadřuje pořadí procesu, na který se result posílá.

Kdybychom použili v předchozím programu příkaz MPI ALLREDUCE místo příkazu MPI REDUCE, pak by se total rozeslal na všechny procesy. Je tedy zřejmé, že syntaxe pří- kazu MPI ALLREDUCE je totožná se syntaxí příkazu MPI REDUCE s tím, že neobsahuje argument root.

2.8 Často používané MPI příkazy

Uveďme si pro přehlednost probrané MPI příkazy a doplňme si je o nové.

2.8.1 Příkazy pro vytvoření a správu paralelního prostředí

A. Inicializace MPI

call MPI_INIT(ierr)

B. Ukončení MPI

call MPI_FINALIZE(ierr)

C. Definování komunikátoru a určení celkového počtu procesů

call MPI_COMM_SIZE(comm,nprocs,ierr)

11

Uvádíme jen nejpoužívanější operace.

(31)

D. Definování pořadí jednotlivých procesů

call MPI_COMM_RANK(comm,myrank,ierr)

E. Zastavení paralelního programu

call MPI_ABORT(comm,errorcode,ierr)

Příkaz MPI ABORT je paralelní verzí příkazu programovacího jazyka STOP. Použije-li jakýkoliv proces příkaz MPI ABORT (např. dojde-li k chybě během čtení na procesu), ukončí se běh paralelního programu na všech procesech.

2.8.2 Příkazy pro kolektivní komunikaci

A. Synchronizace či blokace paralelního programu

call MPI_BARRIER(comm,ierr)

B. Rozeslání informace na jednotlivé procesy

call MPI_BCAST(broadcast_message,count,MPI_data_type,root,comm,ierr)

C. Distribuce informací stejné velikosti na jednotlivé procesy

call MPI_SCATTER(sendbuf,sendcount,send_MPI_data_type, &

recvbuf,recvcount,recv_MPI_data_type,root,comm,ierr)

Funkce příkazu je vysvětlena na obr. 6 a použití MPI SCATTER demonstruje násle- dující program.

program scatter implicit none

!

integer :: i

!

real :: sendmsg(4), & ! send message recvmsg ! recieve message

!

! start up MPI

!

(32)

call MPI_INIT(ierr)

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatter: nprocs > 4!"

!

if(myrank == 0) then do i=1,nprocs

sendmsg(i)=real(i) enddo

endif

!

call MPI_SCATTER(sendmsg,1,MPI_REAL,recvmsg,1,MPI_REAL,0, &

!

print*," myrank =",myrank," recvmsg =",recvmsg

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program scatter

Obrázek 6: Schematický popis funkce příkazu MPI SCATTER.

D. Shromáždění informací stejné velikosti z jednotlivých procesů

call MPI_GATHER(sendbuf,sendcount,send_MPI_data_type, &

recvbuf,recvcount,recv_MPI_data_type,root,comm,ierr)

Funkce příkazu je vysvětlena na obr. 7 a použití MPI GATHER demonstruje následující

program.

(33)

program gather implicit none

!

real :: sendmsg, & ! send message recvmsg(4) ! recieve message

!

! start up MPI

!

call MPI_INIT(ierr)

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gather: nprocs > 4!"

!

sendmsg=real(myrank+1)

!

call MPI_GATHER(sendmsg,1,MPI_REAL,recvmsg,1,MPI_REAL,0, &

!

if(myrank == 0) then print*," recvmsg:",recvmsg endif

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program gather

processes in comm need to call this routine.

Obrázek 7: Schematický popis funkce příkazu MPI GATHER.

E. Distribuce informací nestejné velikosti na jednotlivé procesy

(34)

call MPI_SCATTERV(sendbuf,sendcounts,displs,send_MPI_data_type, &

recvbuf,recvcount, recv_MPI_data_type, &

root,comm,ierr)

Funkce příkazu je vysvětlena na obr. 8 a použití MPI SCATTERV demonstruje násle- dující program.

program scatterv implicit none

!

integer :: send_count(0:3), & ! send count recv_count, & ! receive count displ(0:3) ! displacement

!

real :: sendmsg(10), & ! send message recvmsg(4) ! recieve message

!

! start up MPI

!

call MPI_INIT(ierr)

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatterv: nprocs > 4!"

!

if(myrank == 0) then

sendmsg=(/1.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,4.0,4.0/) send_count=(/1,2,3,4/)

displ=(/0,1,3,6/) endif

recv_count=myrank+1

!

call MPI_SCATTERV(sendmsg,send_count,displ,MPI_REAL, &

recvmsg,recv_count, MPI_REAL, &

0,MPI_COMM_WORLD,ierr)

!

print*," myrank =",myrank," recvmsg:",recvmsg

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program scatterv

(35)

Obrázek 8: Schematický popis funkce příkazu MPI SCATTERV.

F. Shromáždění informací nestejné velikosti z jednotlivých procesů

call MPI_GATHERV(sendbuf,sendcount, send_MPI_data_type, &

recvbuf,recvcounts,displs,recv_MPI_data_type, &

root,comm,ierr)

Funkce příkazu je vysvětlena na obr. 9 a použití MPI GATHERV demonstruje následující program.

program gatherv implicit none

!

integer :: nprocs, & ! # of processes myrank, & ! my process rank ierr,i

integer :: send_count, & ! send count recv_count(0:3), & ! receive count displ(0:3) ! displacement

!

! start up MPI

!

call MPI_INIT(ierr)

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " gatherv: nprocs > 4!"

!

(36)

!

do i=1,myrank+1 sendmsg(i)=myrank+1 enddo

send_count=myrank+1

!

recv_count=(/1,2,3,4/) displ=(/0,1,3,6/)

!

call MPI_GATHERV(sendmsg,send_count, MPI_REAL, &

recvmsg,recv_count,displ,MPI_REAL, &

0,MPI_COMM_WORLD,ierr)

!

if(myrank == 0) then print*," recvmsg:",recvmsg endif

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program gatherv

Obrázek 9: Schematický popis funkce příkazu MPI GATHERV.

Příkazy MPI GATHER a MPI GATHERV mají ještě „ALL” verze MPI ALLGATHER a MPI ALLGATHERV, které shromažďují informace na všechny procesy a tudíž neobsahují argument root.

G. Posílání informací stejné velikosti z jednotlivých procesů na všechny procesy

call MPI_ALLTOALL(sendbuf,sendcount,send_MPI_data_type, &

recvbuf,recvcount,recv_MPI_data_type, &

comm,ierr)

(37)

Funkce příkazu je vysvětlena na obr. 10 a použití MPI ALLTOALL demonstruje násle- dující program.

Obrázek 10: Schematický popis funkce příkazu MPI ALLTOALL.

program alltoall implicit none

!

! start up MPI

!

call MPI_INIT(ierr)

!

do i=1,nprocs

sendmsg(i)=REAL(i+nprocs*myrank) enddo

print*," myrank =",myrank," sendmsg:",sendmsg

!

call MPI_ALLTOALL(sendmsg,1,MPI_REAL, &

recvmsg,1,MPI_REAL, &

!

! shut down MPI

call MPI_FINALIZE(ierr) end program alltoall

(38)

H. Posílání informací nestejné velikosti z jednotlivých procesů na všechny procesy

call MPI_ALLTOALLV(sendbuf,sendcounts,sdispls,send_MPI_data_type, &

recvbuf,recvcounts,rdispls,recv_MPI_data_type, &

comm,ierr)

Funkce příkazu je vysvětlena na obr. 11 a použití MPI ALLTOALLV demonstruje ná- sledující program.

program alltoallv implicit none

!

integer :: scnt(0:3), & ! send counts sdsp(0:3), & ! send displacements rcnt(0:3), & ! recv counts rdsp(0:3) ! recv displacements real :: sendmsg(10),recvmsg(16)

!

sendmsg=(/1.0,2.0,2.0,3.0,3.0,3.0,4.0,4.0,4.0,4.0/) scnt=(/1,2,3,4/)

sdsp=(/0,1,3,6/)

!

! start up MPI

!

call MPI_INIT(ierr)

!

! sendmsg

!

do i=1,10

sendmsg(i)=sendmsg(i)+REAL(nprocs*myrank) enddo

!

! rcnt and rdsp

!

do i=0,nprocs-1 rcnt(i)=myrank+1 rdsp(i)=i*(myrank+1) enddo

!

call MPI_ALLTOALLV(sendmsg,scnt,sdsp,MPI_REAL, &

recvmsg,rcnt,rdsp,MPI_REAL, &

!

(39)

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program alltoallv

Obrázek 11: Schematický popis funkce příkazu MPI ALLTOALLV.

2.8.3 Příkazy pro operace na proměnných distribuovaných na jednotlivých procesech

A. Příkaz „Reduce”

call MPI_REDUCE(operand,result,count,MPI_data_type,operation,root, &

comm,ierr)

Příkaz MPI REDUCE má ještě „ALL” verzi, které neobsahuje argument root.

B. Příkaz „Scan”

call MPI_SCAN(sendbuf,recvbuf,count,MPI_data_type,operation,comm,ierr)

Funkce příkazu je vysvětlena na obr. 12 a použití MPI SCAN demonstruje následující

program.

(40)

program scan implicit none

!

real :: sendmsg,recvmsg

!

! start up MPI

!

call MPI_INIT(ierr)

!

! sendmsg

!

sendmsg=REAL(myrank+1)

!

call MPI_SCAN(sendmsg,recvmsg,1,MPI_REAL,MPI_SUM, &

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program scan

Obrázek 12: Schematický popis funkce příkazu MPI SCAN.

C. Příkaz „Reduce Scatter”

call MPI_MPI_REDUCE_SCATTER(sendbuf,recvbuf,recvcounts,MPI_data_type, &

operation,comm,ierr)

(41)

Funkce příkazu je vysvětlena na obr. 13 a použití MPI REDUCE SCATTER demon- struje následující program.

Obrázek 13: Schematický popis funkce příkazu MPI REDUCE SCATTER.

program reduce_scatter implicit none

!

integer :: rcnt(0:3) ! recv counts

!

rcnt=(/1,2,3,4/)

!

! start up MPI

!

call MPI_INIT(ierr)

!

call MPI_COMM_SIZE(MPI_COMM_WORLD,nprocs,ierr) if(nprocs > 4) stop " scatter: nprocs > 4!"

!

(42)

!

do i=1,10

sendmsg(i)=real(i+myrank*10) enddo

!

call MPI_REDUCE_SCATTER(sendmsg,recvmsg,rcnt,MPI_REAL,MPI_SUM, &

!

print*," myrank =",myrank," recvmsg =",recvmsg

!

! shut down MPI

!

call MPI_FINALIZE(ierr) end program reduce_scatter

(43)

3 APLIKACE

3.1 Paralelní „tempering”

Paralelní „tempering” je způsob jak zefektivnit, urychlit či zlepšit vzorkování např.

v Monte Carlo (MC) metodě.

¹²

Ukážeme si princip paralelního „temperingu” na jednodu- chém případu částice v jednorozměrném silovém poli.

3.1.1 Částice v jednorozměrném silovém poli

Představme si, že máme částici v jednorozměrném silovém poli. Částice má teplotu T a silové pole je charakterizované potenciálem U (x):

U (x) = ∞ x < −2

= 1 · [1 + sin (2πx)] − 2 ≤ x ≤ −1.25

= 2 · [1 + sin (2πx)] − 1.25 ≤ x ≤ −0.25

= 3 · [1 + sin (2πx)] − 0.25 ≤ x ≤ 0.75 (10)

= 4 · [1 + sin (2πx)] 0.75 ≤ x ≤ 1.75

= 5 · [1 + sin (2πx)] 1.75 ≤ x ≤ 2

= ∞ x > 2

Průběh U (x) je graficky zobrazen na obr. 14. Potenciál U (x) je charakterizován čtyřmi lokálními minimy oddělenými bariérami vzrůstající velikosti ve směru kladné osy x.

12

Význam anglického výrazu „tempering” je zřejmý z dalšího výkladu.

(44)

x

-2 -1 0 1 2

U(x)

0 2 4 6 8 10

Obrázek 14: Průběh potenciálu U (x).

3.1.2 Monte Carlo metoda

„Pohyb” částice v intervalu x ∈ h−2; 2i (vzorkování fázového prostoru) budeme reali- zovat pomocí MC metody. MC metoda používá následující algoritmus:

1. Nechť se částice nachází v poloze x

^o

a má zde hodnotu potenciálu U(x

^o

).

2. Vygenerujme novou pozici x

ⁿ

příkazem

x

ⁿ

= x

^o

+ (2ζ − 1) ∆x

_max

(11)

kde ζ je náhodné číslo z intervalu (0; 1) a ∆x

_max

je maximální posunutí ve směru x.

Současně spočtěme hodnotu potenciálu v x

ⁿ

, U(x

ⁿ

).

3. Spočtěme faktor

exp {−β [U (x

ⁿ

) − U (x

^o

)]} (12) a proveďme test zda

exp {−β [U (x

ⁿ

) − U (x

^o

)]} ≥ ζ (13)

V rovnicích (12) a (13) β = 1/(k

B

T ) a k

B

je Boltzmannova konstanta. Pokud je

podmínka (13) splněna, pak nahradíme „starou” pozici x

^o

„novou” pozicí x

ⁿ

. Pokud

podmínka (13) není splněna, „pohyb” se nepřijme a pozice x

^o

se nezmění.

(45)

Výše popsaný postup se opakuje dostatečně dlouho než dojde k náležitému „provzorkování”

fázového prostoru.

MC procesem spočteme pravděpodobnost nalezení částice v místě x, P (x), pro různé teploty T . P (x) určíme následujícím způsobem:

1. Nadefinujeme velikost BIN u pro histogram, ∆x.

2. Příkazem

BIN = INT ((x

^o

− x

_min

)/∆x) + 1 (14) spočteme pořadí BIN u v histogramu. V rovnici (14) je BIN deklarován jako INTE- GER, x

_min

je minimální hodnota x, kterou částice může dosáhnout (v našem případě x

_min

= −2) a INT je vnitřní FORTRANská funkce, jež převádí čísla typu REAL na čísla typu INTEGER.

3. Akumulujeme histogram příkazem

histogram (BIN ) = histogram(BIN ) + 1 (15) 4. Po skončení MC procesu normalizujeme histogram celkovým počtem MC pokusů.

MC program pro částici v jednorozměrném silovém poli může vypadat následovně:

!

! A Single Particle on a 1D Potential with Multiple Local Minima

!

program SerPT implicit none

integer,parameter :: maxhist=1000

real(8),parameter :: x_min=-2.0d0,x_max=2.0d0 integer :: idum,ncyclus,nprint, &

nacp_loc,ntri_loc,nc,ih, &

histogram(maxhist) real(8) :: temp,beta, &

dxmax,x,Ux, &

dx,rat_loc

!

! Auxiliary Variables

!

dxmax=0.1d0 dx=0.025d0

!

! Init idum

!

idum=-471

!

! Input

!

print*," # MC Cyclus:"

read(*,*) ncyclus

print*," Print Frequency:"

read(*,*) nprint

!

! Iput Values (Hard Wired)

!