PFortran w praktyce

Daniel Rychcik
muflon@mat.uni.torun.pl
22-05-2001

Przypomnienie

PFortran jest rozszerzeniem składni Fortranu 77 o dwa operatory służące do komunikacji między procesami:

Procesy w PFortranie numerowane są od 0. Ilość procesów określa zmienna nProc, numer aktualnego procesu - zmienna myProc. Dostępna jest również zmienna nProc0 równa nProc-1.

PFortran wprowadza również kilka dodatkowych procedur z których dla użytkownika dostępne są:

Translator języka działa w ten sposób, że kod w PFortranie (plik z rozszerzeniem .pf) jest konwertowany do postaci zwykłego programu w Fortranie 77 używającego wywołań biblioteki MPI. Następnie program jest kompilowany kompilatorem f77 i linkowana jest z nim biblioteka pfmpi (oraz odpowiednia biblioteka MPI dla Fortranu).

Translator napisany jest w języku C z użyciem generatora analizatorów składniowych bison. Implementacja jego rozszerzeń będzie prawdopodobnie w dużej części polegała na rozszerzeniu plików źródłowych bisona.

Cechy

Główną zaletą PFortranu w porównaniu z bibliotekami równoległymi typu MPI czy PVM jest prostota i zwięzłość kodu. Pozwala to na pisanie programów w sposób bardziej czytelny - publikacje dotyczące algorytmów równoległych używają podobnej notacji do ich zapisu. Daje to również szansę na pisanie programów równoległych ludziom, którzy nie mają ochoty uczyć się skomplikowanej składni wywołań funkcji bibliotek message passing mając do napisania naprawdę prosty program.

Kod w PFortranie jest przenośny - wszystkie szczegóły implementacyjne są ukryte w bibliotece dołączanej przy linkowaniu oraz w generowanym kodzie Fortranu 77.

Niestety uproszczenie modelu programowania pociąga za sobą pewien spadek wydajności w porównaniu z ręcznym przesyłaniem danych - w tym momencie nie mamy już wpływu na to, jak dane będą przesyłane. Dalej podaję przykład takiej sytuacji (operacja redukcji).

Szczegóły techniczne

Zaawansowana redukcja, własne operatory

Operacja redukcji w PFortranie może zostać zapisana na dwa sposoby: jako operacja przypisania:

       wynik=FUN{op}

lub jako wywołanie procedury dwuargumentowej:

       call FUN{wynik,op}

W drugim przypadku wynik redukcji zostanie zwrócony w pierwszym argumencie jako efekt uboczny działania procedury.

Funkcja będąca argumentem operacji redukcji musi być przemienna. Może nią być któraś z wbudowanych funkcji Fortranu (np. AMIN), lub własna funkcja użytkownika (patrz niżej). Kod jest generowany w taki sposób, że na danej maszynie przy takiej samej ilości procesorów wyniki będą powtarzalne (kolejność wykonywania działań jest funkcją ilości procesorów).

Własne operatory redukcji możemy tworzyć definiując funkcje lub procedury o następującej postaci:

  <typ> FUNCTION <nazwa>(<arg1>,<arg2> [,<dodatkowe argumenty>]) 
  SUBROUTINE <nazwa>(<arg1>,<arg2>,<wynik> [,<dodatkowe argumenty>])

W przypadku definiowania procedury musi ona zwracać wynik w trzecim argumencie. Funkcje redukcji powinny być typu REAL lub INTEGER. Dodatkowe parametry funkcji/procedury są opcjonalne i przy wywołaniu będą traktowane jako globalne - identyczne na wszystkich procesorach.

Parametry procedury (arg1,arg2 i wynik) muszą być tej samej postaci - jeżeli dane są tablicami dwuelementowymi, to wynik musi być również tablicą dwuelementową.

Oto przykładowy prosty operator redukcji (własna suma):

       REAL*8 FUNCTION MYSUM(IN1,IN2)
         MYSUM=IN1+IN2
         RETURN
       END

Możemy go używać podobnie jak standardowego: WYNIK=MYSUM{ELEM}.

Sposób definicji operatora redukcji (procedura/funkcja) nie ma nic wspólnego ze sposobem jego używania! Oznacza to, że możemy np. zdefiniować sumowanie kilkuelementowych tablic jako procedurę, a używać jej następnie przy pomocy składni TAB=SUMA{TAB}. Upraszcza to trochę programowanie, bo nie musimy definiować tablic pomocniczych do przekazywania danych.

Wbudowany operator redukcji +{...} pozwala na podawanie jako argumentów tablic jednowymiarowych. Przykładowo następujący kod wyliczy sumy elementów na każdej pozycji tablicy A i umieści te sumy w tablicy B na wszystkich procesorach:

       B=+{A}

"Message tag"

Ograniczenia translatora PFortranu (opisane niżej) są efektem sposobu realizacji komunikacji wybranego przez autorów. Każdy komunikat jest bowiem opatrywany unikalnym znacznikiem (tag). Jest on zwiększany o 1 po każdej operacji komunikacji i używany przez procedury komunikacyjne niskiego poziomu (np. MPI) do rozróżniania komunikatów. Jeżeli zatem na niektórych procesach nie wykona się instrukcja postaci a=b@1 to liczniki znaczników się rozsynchronizują i dalsza komunikacja będzie niemożliwa!

Po obejrzeniu kodu generowanego np. przez operator redukcji widać, że zwiększenie licznika realizuje procedura pf_inctag()

Operacje wejścia/wyjścia

I/O w PFortranie można realizować na kilka sposobów, przy czym wybór konkretnej realizacji jest uzależniony z reguły od warunków sprzętowych którymi dysponujemy.

Pierwsza, najczęściej stosowana metoda jest taka, że wyróżniony proces (najczęściej proces numer 0) czyta stosowne dane z dysku i wysyła je do pozostałych procesów używająć standardowych sposobów komunikacji PFortranu. Oto przykład:

       IF (myProc.EQ.0) CALL OpenFile(coord,unit)
       DO 100,P=0,nProc-1
         IF (myProc.EQ.0) READ(unit,*)(X(i),i=1,N/P)
         X(1:N/P)@p=X(1:N/P)@(nProc-1)
 100   CONTINUE
       IF (myProc.EQ.0) CLOSE(unit)

Inną metodą jest wczytywanie danych kolejno bezpośrednio przez zainteresowane procesy. Należy wtedy zadbać o ich synchronizację, aby nie powodować zbędnego przeciążenia systemu dyskowego:

       DO 110,P=0,nProc-1
         IF (P.EQ.myProc) THEN
           CALL OpenFile(coord,unit)
           DO 100,B=0,myProc-1
             READ(unit,*)(X(I),I=1,N/P)
 100       CONTINUE
           READ(unit,*)(X(i),i=1,N/P)
           CLOSE(unit)
         ENDIF
 110   CONTINUE
       CALL Synchronize

Jeżeli dysponujemy odpowiednim sprzętem, umożliwiającym równoległe wczytywanie danych z pliku do różnych procesów, możemy pokusić się o następujące rozwiązanie:

       CALL OpenFile(coord,unit)
       CALL PositionFile(unit,(N/P)*myProc-1)
       READ(unit,*)(X(I),I=1,N/P)
       CLOSE(unit)

Pomocnicze operacje wejścia/wyjścia takie jak wydruki diagnostyczne, możemy wykonywać na podobnych zasadach, jak z użyciem samej biblioteki message-passing. W przypadku MPICH oznacza to, że komunikaty pojawią się na ekranie w kolejności numerów procesów, które je wysłały.

Ograniczenia translatora

1. Instrukcje zawierające operatory @ i {} muszą się zawsze wykonywać równocześnie na wszystkich procesorach

Nieprawidłowy jest zatem kod zawierający odwołania do nich w instrukcji warunkowej której warunek jest uzależniony od myProc:

       IF (myPROC.EQ.0) THEN
         suma=+{suma}
       ENDIF

Nie jest to jednak dużym ograniczeniem, gdyż to samo możemy zapisać:

       suma@0=+{suma}

(Ale patrz ograniczenie numer 3!)

2. Instrukcje przypisania zawierające @ muszą mieć na każdym procesorze identyczne adresy źródłowe i docelowe

Należy zatem szczególnie uważać, jeśli adresy te nie są stałe - nie możemy na przykład napisać:

       p=0
       IF (myProc.EQ.0) THEN
         p=1
       ENDIF
       j=j@p

ponieważ proces zerowy będzie miał w ostatniej linii inną wartość p. Zamiast tego możemy napisać na przykład tak:

       j=j@0
       j@0=j@1

3. Złożone wyrażenia zawierające @ i {}

W obecnej wersji translator nie radzi sobie najlepiej ze złożonymi wyrażeniami zawierającymi nowe operatory. Należy zatem starać się pisać kod prosty, w którym przesyłanie danych jest realizowane etapami.

4. Użycie operatora @

Jak na razie, użycie @ jest ograniczone do instrukcji przypisania. Nie powinien on zatem pojawiać się np. w warunkach pętli, lub argumentach funkcji.

5. Specjalne zasady używania myProc

Jeśli chcemy używać zmiennej myProc w adresach operatora @, to musimy ją wpisać bezpośrednio, bez używania zmiennych pomocniczych. Zatem zamiast pisać:

       p=myProc
       a@p=b@(p-1)

powinniśmy pisać:

       a@myProc=b@(myProc-1)

pamiętając oczywiście o zastrzeżeniu numer 1

Zauważone mankamenty

Rozsyłanie tablic wielowymiarowych

Pisząc program na równoległe mnożenie macierzy (matmul.pf) zauważyłem, że translator generuje nie do końca działający kod w przypadku używania operatora @ do przesyłania większych fragmentów tablic dwuwymiarowych. Oto przykładowy kod:

         DO 302,I=0,nProc-1
           DO 301,J=PROC_START(I),PROC_START(I+1)-1
             C(J,1:N)=C(J,1:N)@I 
 301       CONTINUE
 302     CONTINUE

W założeniu ten fragment miał rozsyłać obliczone fragmenty macierzy wynikowej do wszystkich procesów. Niestety po przetłumaczeniu kodu na zwykły FORTRAN wynik translacji instrukcji C(J,1:N)=C(J,1:N)@I wygląda następująco:

      o10000=((n-1+1))*4
      o10001=1
      do 55000 o10002 = 1, n
        pf_rb(o10001)=c(j,o10002)
        o10001=o10001+1
55000 continue
        call pf_brdst(i, pf_rb(1), o10000)

Widać tu mniej więcej intencję translatora: budujemy wektor pf_rb z odpowiednich elementów macierzy C i wykonujemy broadcast w taki sposób, że wszystkie procesy otrzymają fragment zbudowany przez proces I. Niestety brakuje tu jeszcze jednej instrukcji - przepisania otrzymanego wektora pf_rb do odpowiedniego fragmentu macierzy C! Po ręcznym wstawieniu odpowiedniej pętli program oczywiście działa, jednak z pewnością nie jest to dobre rozwiązanie.

Widać też, że translator nie potrafi skorzystać ze znajomości rozmieszczenia elementów macierzy w tablicy - tak naprawdę, to pomocniczy wektor nie jest tu wcale potrzebny! Aby ten kod jednak zadziałał, należy go w PFortranie zapisać następująco:

         DO 303,I=0,nProc-1
           DO 302,J=PROC_START(I),PROC_START(I+1)-1
             DO 300,K=1,N
               ROW(K)=C(J,K)
 300         CONTINUE
             ROW(1:N)=ROW(1:N)@I
             DO 301,K=1,N
               C(J,K)=ROW(K) 
 301         CONTINUE
 302       CONTINUE
 303     CONTINUE

Ale to jest wlasciwie przepisywanie kodu generowanego przez translator!

Implementacja redukcji

Po obejrzeniu wyniku kompilacji programu dotprod.pf widać, że do wyliczenia sumy stosowany jest algorytm drzewa binarnego, jednak jest ona wyliczana najpierw w procesie 0, a następnie rozsyłana do pozostałych procesów (również z użyciem struktury drzewiastej). O wiele efektywniejszy byłby tutaj algorytm hypercube - jest to przeciez odpowiednik operacji AllReduce standardu MPI. Oto kod generowany dla prostej instrukcji SUM=+{SUM}:

Cpf          SUM=+{SUM}
      sum=sum
      sum=sum
Cpf
Cpf   dimxxx: non-powers of 2 allowed
      o10001=pflog2ceil(nProc)
Cpf
Cpf   reduce to process 0
      do 55000 o10006 = o10001, 1,-1
        o10005=2**o10006
        o10004=2**(o10006-1)
        call pf_inctag()
        if (myProc.LT.o10004) then
          o10003=xor(myProc,o10004)
          if (o10003.LE.nProc0) then
            call pf_rcv(o10000,4)
            sum=sum+o10000
          endif
        else if (myProc.LT.o10005) then
          o10002=xor(myProc,o10004)
          call pf_snd(o10002,sum,4)
        endif
55000 continue
Cpf
Cpf   propogate result up from zero
      do 55001 o10006 = 1, o10001
        o10004=2**o10006
        o10005=2**(o10006-1)
        call pf_inctag()
        if (myProc.LT.o10005) then
          o10002=xor(myProc,o10005)
          if (o10002.LT.nProc) then
            call pf_snd(o10002,sum,4)
          endif
        else if (myProc.LT.o10004) then
          o10003=xor(myProc,o10005)
          call pf_rcv(sum,4)
        endif
55001 continue

Co dziwne, zauważyłem drobne rozbieżności w wynikach w zależności od ilości procesów! Nie można tego zrzucić na błędy zaokrągleń, gdyż pojawiały się one nawet przy wektorach złożonych z kilkudziesięciu elementów typu REAL*8.



Literatura

  1. B.Bagheri, T.W.Clark, L.Ridgway Scott: PFortran Reference Manual
    http://planguages.cs.uchicago.edu/docs/pfmanual.ps.
  2. S.Petrykowski: PFortran
    http://www.mat.uni.torun.pl/~bala/sem_mgr_2000/pf_axon.html.