Entwicklung eines parallelen Rahmenprogramms

Eine Parallelisierung bezüglich der Komponenten des gewöhnlichen Differentialgleichungssystems kann für jede der in Abschnitt2.4vorgestellten Implementierungsvarianten und auch für viele alternative Rea-lisierungen des Berechnungskernel durchgeführt werden, indem die Iterationen der Schleifen über die Sy-stemdimension parallel durch mehrere Kontrollflüsse ausgeführt werden. Dazu ist es zunächst erforderlich, eine geeignete Verteilung der Komponenten bzw. der zu den Komponenten korrespondierenden Schleifen-iterationen festzulegen. Zusätzlich müssen geeignete Synchronisationsoperationen eingefügt werden, die sicherstellen, daß die aus der parallelen Ausführung hervorgehende Berechnungsreihenfolge in Überein-stimmung mit den durch die Berechnungsvorschrift (2.2) gegebenen Datenabhängigkeiten steht.

3.2.1 Wahl der Datenverteilung

Da die Komponentenfunktionenf_junabhängig voneinander ausgeführt werden können und deshalb keine der Schleifen über die Systemdimension Abhängigkeiten zwischen den Iterationen der Schleifen aufweist, ist prinzipiell eine beliebige Verteilung der Komponenten auf die beteiligten Kontrollflüsse möglich. Die Wahl der Datenverteilung kann sich jedoch auf das Lokalitätsverhalten und den zur Synchronisation der Kontrollflüsse erforderlichen Aufwand auswirken. In Frage kommen z. B. typische Datenverteilungen, wie sie in vielen Anwendungen aus dem wissenschaftlichen Rechnen für die Verteilung eindimensionale Felder verwendet werden (vgl.Rauber u. Rünger 2000). Alle diese Unterteilungen haben zum Ziel, jedem dernK

Kontrollflüsse etwa den gleichen Arbeitsaufwand zuzuordnen, indem die Indexmenge der Komponenten I={j∈N: 1≤j≤n}innKannähernd gleich große TeilmengenI_k,k=1, . . . ,nK, partitioniert wird.

3.2 Entwicklung eines parallelen Rahmenprogramms 93 I Blockweise Verteilung:Die MengeIwird durch Festlegung vonn_K−1 Intervallgrenzen

j₁<j₂<· · ·<j_nK−1 (3.1) innKdisjunkte Intervalle (Blöcke)

I_k = (jk−1,j_k] ={j∈I:jk−1<j≤ j_k} für k=1, . . . ,nK (3.2) unterteilt, wobeij0 = 0 und jn_K = n gesetzt wird. Zur Festlegung der Intervallgrenzen kommen ver-schiedene Vorgehensweisen in Frage. Der einfachste Fall liegt vor, wennnein ganzzahliges Vielfaches von nK ist. Dann können jedem Kontrollfluß genau n/nK Komponenten zugeordnet werden, indem j_k = k·n/nKgesetzt wird. Istnkein ganzzahliges Vielfaches vonnK, legt man die Intervallgrenzen so fest, daß(nmodnK)Kontrollflüssedn/nKeKomponenten erhalten und allen verbleibenden Kontroll-flüssenbn/nKcKomponenten zugeordnet werden, z. B.:

j_k=

I Zyklische Verteilung:Die Elemente der IndexmengeIwerden reihum an die beteiligten Kontrollflüsse verteilt, z. B.:

I_k={j∈ I:(j−1)modnK=k−1} für k=1, . . . ,nK . (3.4) I Blockzyklische Verteilung:Ähnlich einer blockweisen Verteilung wird die Indexmenge I zunächst in

n_Bn_KIntervalle (Blöcke)

B_b= (jb−1,j_b] ={j∈I:jb−1<j≤j_b} für b=1, . . . ,n_B (3.5) mit den Intervallgrenzen

0=j₀<j₁< j₂<· · ·<j_nB−1<jn_B =n (3.6) zerlegt. Diese Blöcke werden anschließend zyklisch den beteiligten Kontrollflüssen zugeordnet:

I_k= ^[

b=1,...,n_B (b−1)modn_K=k−1

B_b für k=_{1, . . . ,}n_K . (3.7)

Um eine möglichst hohe Lokalität der Speicherzugriffe bei der parallelen Implementierung eines einge-betteten Runge-Kutta-Verfahrens zu erreichen, ist in der Regel eine blockweise Datenverteilung die beste Wahl. Da jedem Kontrollfluß ein einzelner zusammenhängender Block von Komponenten zugeordnet wird, ist in diesem Fall die bestmögliche Ausnutzung von zeitlicher und räumlicher Lokalität zwischen den Itera-tionen der Schleifen über die Systemdimension gegeben, wenn der Block komponentenweise durchlaufen wird. Auch für eine effiziente Realisierung des Datenaustauschs ist eine blockweise Verteilung besonders günstig. Eine ausführlichere Diskussion der Vorteile einer blockweisen Verteilung wird später im Zusam-menhang mit der Erstellung von Implementierung für spezifische Speicherarchitekturen durchgeführt.

3.2.2 Synchronisation der Schrittweitenkontrolle

Während des Integrationsvorgangs, d. h. während des Durchlaufens des Integrationsintervalls müssen die Zeitschritte von den beteiligten Kontrollflüssen synchron ausgeführt werden. Das bedeutet, zu jedem Zeit-punkt berechnen alle Kontrollflüsse jeweils den gleichen Zeitschritt und verwenden dazu eine einheitliche Schrittweite. Insbesondere muß eine einheitliche Entscheidung darüber getroffen werden, ob ein Zeitschritt akzeptiert oder verworfen werden soll.

Zur Durchführung der Schrittweitenkontrolle muß zunächst die Approximation des lokalen Fehlers e= ¹

berechnet werden. Da jedem Kontrollflußkaufgrund der Verteilung der Vektorkomponenten nur die durch die IndexmengeI_kbestimmten Komponenten der Vektoreneundsbekannt sind, bietet es sich an, daß alle Kontrollflüsse parallel jeweils ihr lokales Maximume_local^k bestimmen:

e^k_local= ¹ TOLmax

j∈I_k

e[j] s[j]

. (3.9)

Anschließend kann das globale Maximum aller Kontrollflüsse e= max

k=1,...,n_Ke^k_local (3.10)

bestimmt werden, wozu allerdings ein Datenaustausch zwischen den Kontrollflüssen erforderlich ist. Zwei mögliche Vorgehensweisen kommen in Betracht:

1. Nur ein einzelner Kontrollfluß kennt das globale Maximume. Dieser Kontrollfluß führt die Schrittwei-tenkontrolle allein durch, d. h., er entscheidet allein, ob der Schritt akzeptiert oder verworfen wird, und bestimmt die neue Schrittweite. Das Ergebnis der Entscheidung sowie die neue Schrittweite wird an-schließend den anderen Kontrollflüssen bekannt gegeben, die in der Zwischenzeit auf das Eintreffen dieser Daten gewartet haben.

2. Alle Kontrollflüsse kennen das globale Maximume. Anstatt zu warten, können die Kontrollflüsse jetzt die Schrittweitenkontrolle gleichzeitig durchführen. Da alle Kontrollflüsse vom gleichen Wert füre aus-gehen, treffen sie die gleiche Entscheidung bezüglich des Akzeptierens oder Verwerfens des Zeitschrittes und berechnen einen identischen Wert für die neue Schrittweite.

Zur Realisierung der ersten Vorgehensweise müssen an zwei Programmstellen Daten ausgetauscht wer-den. Eine Einzelakkumulationsoperation(auch:Reduktionsoperation) wird benötigt, um anhand der von den einzelnen Kontrollflüssen berechneten lokalen Maximae_local^k das globale Maximumezu bestimmen und dem ausgewählten Kontrollfluß zu übermitteln. Später wird die neue Schrittweite sowie die Information, ob der Zeitschritt akzeptiert oder verworfen wird, mit Hilfe einer Broadcastoperation an die übrigen Kon-trollflüsse weitergegeben. Die zweite Vorgehensweise erfordert lediglich an einer Programmstelle einen Datenaustausch. Hier wird mittels einer spezialisierten Variante einerMultiakkumulationsoperationdas glo-bale Maximumebestimmt und allen Kontrollflüssen bekannt gegeben.

Welche Vorgehensweise effizienter ist, hängt im wesentlichen von der Realisierung des Datenaustauschs ab. Welche Möglichkeiten es dazu gibt, ist wiederum abhängig von der Speicherarchitektur und dem Auf-bau des Verbindungsnetzwerks. Wird die benötigte Funktionalität der Multiakkumulationsoperation bei-spielsweise durch eine Einzelakkumulationsoperation und eine Broadcastoperation implementiert, werden sich die resultierenden Laufzeiten beider Vorgehensweisen in der Regel nur unwesentlich unterscheiden.

Alle in dieser Arbeit betrachteten Implementierungen nutzen die zweite Vorgehensweise.

3.2.3 Realisierung des Rahmenprogramms

Aus den zuvor beschriebenen Überlegungen ergibt sich das in Abb.3.1gezeigte Rahmenprogramm. Das Rahmenprogramm verwendet das sogenannte SPMD-Programmiermodell (englisch:single program multiple data), d. h., alle beteiligten Kontrollflüsse führen das gleiche Programm aus. Anhand der Eingabedaten, beispielsweise der ID des Kontrollflusses, können jedoch gegebenenfalls unterschiedliche Verzweigungen ausgeführt werden, was für die Realisierung des Rahmenprogramms jedoch nicht erforderlich ist.

Es wird eine blockweise Datenverteilung verwendet, die in Zeile 3 des Rahmenprogramms berechnet wird. Dazu bestimmt jeder Kontrollflußkden Index der ersten Komponentej_first=jk−1+1 und den Index der letzten Komponente j_last= j_kseines Datenbereichs I_k = (jk−1,j_k]nach der Vorschrift (3.3). Nach Fest-legung der Datenverteilung wird die Integration vorbereitet, d. h., benötigter Speicher wird allokiert und die Anfangsschrittweite wird bestimmt. Danach wird das Integrationsintervall durch Ausführung der Zeit-schritte durchlaufen. Abschließend gibt der Integrator den nicht mehr benötigten Speicher frei und kehrt zum aufrufenden Programm zurück. Welche Datenstrukturen verwendet und wie sie realisiert werden, hängt von der konkreten Implementierung und von der Speicherarchitektur der Zielplattform ab.