Das Paradigma der parallelen Programmierung und das Rechnen auf Supercomputern
Einführung
In diesem Unterpunkt werden wir uns mit dem Paradigma der parallelen Programmierung befassen. In diesem Programmierparadigma wird bei der Konzeption von rechenintensiven Computerprogrammen berücksichtigt, dass die Rechenleistung des Computers stets voll ausgelastet ist und separate, voneinander unabhängige Teilaufgaben (Tasks) innerhalb der Programme möglichst gleichzeitig (parallel) berechnet werden.
- "Shared Memory" Mehrprozessorsysteme:
Mehrere Prozessoren eines Computers teilen sich einen gemeinsamen Arbeitsspeicher (RAM). Parallelisierung mittels OpenMP (Open Multi-Processing) relativ einfach möglich. - "Distributed Memory" Mehrprozessorsysteme:
Bei Hochleistungs-Großrechenanlagen sind mehrere Mehrprozessor-Computer zu Computer-Clustern verbunden, wobei jeder Computer seinen eigenen privaten Arbeitsspeicher hat. Parallelisierung mittels MPI (Message Passing Interface).
Die Art und Weise wie ein paralleles Programm zu konzipieren ist, hängt auch von der zugrundeliegenden Rechnerarchitektur des ausführenden Computers ab und man unterscheidet hier grob in die nebenstehenden Varianten. Viele Programmiersprachen lassen sich relativ einfach mittels OpenMP (Open Multi-Processing) parallelisieren, welches sich jedoch nur auf "Shared Memory" Systeme anwenden lässt. MPI (Message Passing Interface) stellt dagegen eine weit verbreitete Parallelisierungsmöglichkeit dar, die sowohl auf "Shared Memory" als auch auf "Distributed Memory" Systemen anwendbar ist. Die Parallelisierung eines sequentiellen Programms mittels MPI ist jedoch ein wenig aufwendiger. Dieser Unterpunkt gibt eine Einführung in die parallele Programmierung mit C++ und OpenMP. Die Parallelisierung eines Python-Programms ist aber natürlich auch möglich und kann z.B. mittels des Python-Moduls (threading oder multiprocessing) implementiert werden, bzw. unter Zuhilfenahme von MPI for Python realisiert werden.
Parallelen Programmierung mit OpenMP
Die parallele Erweiterung eines sequentiellen C++ Programms (ein Programm ohne Parallelisierung) ist mittels OpenMP relativ einfach möglich. Betrachten wir z.B. das folgende C++ Programm, welches für neun unterschiedliche q-Werte den Wert der unendlich geometrischen Reihe berechnet (siehe Übungsblatt Nr. 3, Aufgabe 1). Um die Berechnung möglichst genau zu machen, wurde hierbei ein sehr hoher Wert (50 Millionen) als Obergrenze der Partialsumme benutzt \[ \begin{equation} \sum_{k=0}^{50000000} q^k \, \approx \, \sum_{k=0}^\infty q^k = \frac{1}{1-q},\quad \hbox{mit:} \quad q \in ℝ \,\, , \,\, \left| q \right| < 1 \quad . \end{equation} \]
/* Konvergenzverhalten der geometrischen Reihe für unterschiedliche q * Sequentielle Version der Berechnung von 9 approximierten Werten */ #include <iostream> // Ein- und Ausgabebibliothek #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) using namespace std; // Benutze den Namensraum std int main(){ // Hauptfunktion int startTime = time(NULL); // Starten der Zeitmessung double q; // Deklaration des Parameters q const unsigned long int N = 50000000; // Definition der Anzahl der mitgenommenen Summenglieder double wert_approx; // Deklaration des approximierten Wertes der Summe for(int i=1; i <= 9; ++i){ // Schleife zur ueber neun q-Werte q = i*0.1; // Festlegung des q-Wertes wert_approx = 0; // Summenvariable auf Null setzen for(int k=0; k<=N; ++k){ // Schleife zur Berechnung der Summe wert_approx = wert_approx + pow(q,k); // Innerer Teil der geometrischen Reihe } // Ende der Schleife der Summenberechnung printf("q=%6.4f , Wert der Summe = ",q); // Terminalausgabe der berechneten Werte printf("%13.10f (Fehler zum wirklichen Wert : %13.6e ) \n",wert_approx, wert_approx - 1.0/(1-q)); } // Ende der for-Schleife ueber die unterschiedlichen q-Werte // Terminalausgabe der benoetigten Zeit cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; } // Ende Hauptfunktion
Das Programm berechnet nacheinander, in sequentieller Weise die Werte der geometrischen Reihe für neun unterschiedliche q-Werte und gibt die berechneten Werte im Terminal aus. Zusätzlich wurde eine Zeitmessung der benötigten Rechenzeit im Programm implementiert, die am Ende des Programms ebenfalls ausgegeben wird. Die nebenstehende Abbildung zeigt die Terminalausgabe des Programms. Die berechneten Werte stimmen sehr gut mit den wirklichen Werten der unendlich geometrischen Reihe überein, wobei das Programm jedoch zur Berechnung der neun Werte 28 Sekunden benötigt.
Die parallele Erweiterung dieser sequentiellen Version des Programms ist mittels OpenMP relativ einfach möglich. Da die Berechnung der einzelnen Werte nicht voneinander abhängt und somit die neun unterschiedlichen Berechnungen neun unabhängige 'Tasks' darstellen, kann man die for-Schleife, die die Berechnung der einzelnen Werte durchführt, parallel ausführen. Das folgende C++ Programm stellt diese parallele Version dar:
/* Konvergenzverhalten der geometrischen Reihe für unterschiedliche q * Parallele Version der Berechnung von 9 approximierten Werten * Parallelisierung in der for-Schleife ueber neun q-Werte */ #include <iostream> // Ein- und Ausgabebibliothek #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) #include <omp.h> // OpenMP zum parallelen Rechnen using namespace std; // Benutze den Namensraum std int main(){ // Hauptfunktion int startTime = time(NULL); // Starten der Zeitmessung double q; // Deklaration des Parameters q const unsigned long int N = 50000000; // Definition der Anzahl der mitgenommenen Summenglieder double wert_approx; // Deklaration des approximierten Wertes der Summe #pragma omp parallel for private(q,wert_approx) for(int i=1; i <= 9; ++i){ // Schleife zur ueber neun q-Werte q = i*0.1; // Festlegung des q-Wertes wert_approx = 0; // Summenvariable auf Null setzen for(int k=0; k<=N; ++k){ // Schleife zur Berechnung der Summe wert_approx = wert_approx + pow(q,k); // Innerer Teil der geometrischen Reihe } // Ende der Schleife der Summenberechnung printf("q=%6.4f , Wert der Summe = ",q); // Terminalausgabe der berechneten Werte printf("%13.10f (Fehler zum wirklichen Wert : %13.6e ) \n",wert_approx, wert_approx - 1.0/(1-q)); } // Ende der for-Schleife ueber die unterschiedlichen q-Werte // Terminalausgabe der benoetigten Zeit cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; } // Ende Hauptfunktion
Um das Programm zu parallelisieren, bindet man zunächst die OpenMP-Header-Datei mittels '#include <omp.h>' ein. Dann schreibt man vor die zu parallelisierende for-Schleife die Zeile '#pragma omp parallel for '. Dies hat den Effekt, dass nun die einzelnen Summenberechnungen mit den unterschiedlichen q-Werten nicht mehr sequentiell nacheinander ausgeführt werden, sondern mehrere sogenannte Threads die Berechnungen simultan, gleichzeitig ausführen. Spezifiziert man die Anzahl der Threads nicht explizit (mittels omp_set_num_threads()), so entscheidet OpenMP selbständig, in Abhängigkeit von der auf dem Computer zu Verfügung stehenden Prozessor-Kernen, wie viele Threads erzeugt werden sollen. Bei der hier vorliegenden for-Schleife ist bei der Parallelisierung noch das folgende zu beachten: Da die Variablen 'q' und 'wert_approx' von allen Threads benutzt werden, ist es erforderlich diese Variablen beim Parallelisierungsmechanismus als 'private' Thread-Variablen zu kennzeichnen (private(q,wert_approx)), da sie sonst die berechneten Größen der anderen Threads überschreiben. Bemerkung: Man hätte dies auch vermeiden können, falls man die beiden Variablen lokal innerhalb der for-Schleife deklariert hätte.
Um ein OpenMP-Programm auch parallel ausführen zu können, muss man es mittels des Zusatzes '-fopenmp' kompilieren ( g++ -fopenmp Geom_Reihe_parallel.cpp ). Die nebenstehende Terminalausgabe zeigt, dass das Programm (auf meinem Laptop) nur 5 Sekunden benötigt. Des Weiteren erkennt man, dass sich die Reihenfolge der ausgegebenen Werte, im Vergleich zum sequentiellen Programm, geändert hat. Die Ausgabe erfolgt nun nicht mehr in Reihenfolge, sondern jeder Thread gibt nach Beendigung seines Tasks den berechneten Wert im Terminal aus. Die nebenstehende Ausgabe zeigt die Auslastung der einzelnen Prozessorkerne (CPUs) meines Laptops beim Ausführen der sequentiellen und parallelen Version des Programms. Beim Ausführen des sequentiellen Programms wird nur ein Prozessorkern benutzt, wohingegen bei der parallelen Version viele der verfügbaren CPUs simultan die einzelnen Werte berechnen und somit viel weniger Zeit benötigt wird (5 Sekunden anstatt 28 Sekunden).
Man hätte die Parallelisierung des Programms auch auf die for-Schleife der Summenberechnung anwenden können, wie die folgende Version des parallelen Programms zeigt.
/* Konvergenzverhalten der geometrischen Reihe für unterschiedliche q * Parallele Version der Berechnung von 9 approximierten Werten * Parallelisierung in der for-Schleife zur Berechnung der Summe bei festem q */ #include <iostream> // Ein- und Ausgabebibliothek #include <cmath> // Bibliothek für mathematisches (e-Funktion, Betrag, ...) #include <omp.h> // OpenMP zum parallelen Rechnen using namespace std; // Benutze den Namensraum std int main(){ // Hauptfunktion int startTime = time(NULL); // Starten der Zeitmessung double q; // Deklaration des Parameters q const unsigned long int N = 50000000; // Definition der Anzahl der mitgenommenen Summenglieder double wert_approx; // Deklaration des approximierten Wertes der Summe for(int i=1; i <= 9; ++i){ // Schleife zur ueber neun q-Werte q = i*0.1; // Festlegung des q-Wertes wert_approx = 0; // Summenvariable auf Null setzen #pragma omp parallel for reduction(+:wert_approx) for(int k=0; k<=N; ++k){ // Schleife zur Berechnung der Summe wert_approx = wert_approx + pow(q,k); // Innerer Teil der geometrischen Reihe } // Ende der Schleife der Summenberechnung printf("q=%6.4f , Wert der Summe = ",q); // Terminalausgabe der berechneten Werte printf("%13.10f (Fehler zum wirklichen Wert : %13.6e ) \n",wert_approx, wert_approx - 1.0/(1-q)); } // Ende der for-Schleife ueber die unterschiedlichen q-Werte // Terminalausgabe der benoetigten Zeit cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; } // Ende Hauptfunktion
Das Programm benötigt ebenfalls nur 5 Sekunden um die neun Werte zu berechnen und hat zusätzlich den Vorteil, dass die Terminalausgabe in gewohnter Reihenfolge ausgegeben wird. Bei der Parallelisierung ist nun eine Kennzeichnung von privaten Variablen nicht mehr nötig, jedoch ist der Zusatz 'reduction(+:wert_approx)' erforderlich. In der hier gewählten Parallelisierungsvariante werden die einzelnen Teile der Summenberechnung auf mehrere Threads aufgeteilt und es kann dabei ein Schreib-Lese-Konflikt (‘race condition’) in der Zeile 'wert_approx = wert_approx + pow(q,k);' auftreten; der Zusatz 'reduction(+:wert_approx)' verhindert diesen Konflikt.
Beispiel: Die parallele OpenMP-Version des C++ Programms "Eine Kiste voller Pendel"
Im Folgenden werden wir eine parallele OpenMP-Version des C++ Programms "Eine Kiste voller Pendel" (siehe Vorlesung 11, Unterpunkt Klausurvorbereitung und prüfungsrelevante Themen) erstellen. Wir benutzen dafür die folgende sequentielle Version des C++ Programms, wobei wir im Vergleich zur ursprünglichen Version zusätzlich eine Zeitmessung eingebaut haben und die Anzahl der Zeitgitterpunkte und das zu simulierende Zeitintervall erhöht haben:
/* Sequenzielle Version * Beispiel einer, von der Basisklasse 'Ding' abgeleitete Klasse 'Pendel' * Zwei zusaetzliche private Instanzvariablen (l und beta) und die winkelspezifischen Anfangswerte (phi und v_phi) * wurden der abgeleiteten Pendel-Klasse hinzugefuegt * Der Konstruktor der Klasse Pendel erzeugt ein Objekt Ding und initialisiert seine eigenen Daten-Member * Von der Klasse Pendel wurde eine weitere Subklasse 'Pendel_math' abgeleitet, welche die lineare Näherung der Pendel-DGL benutzt (gedaempfter harmonischer Ozillator) * Direkte Ausgabe der berechneten Werte von 25 Pendel-Simulationen in eine Datei (Visualisierung mittels Pendel_Container.py) * Zusatzliche Zeitmessung zum späteren Vergleich mit der parallelen Variante */ #include "Pendel.hpp" // Pendel und Ding Klassen #include <iostream> // Ein- und Ausgabebibliothek #include <vector> // Sequenzieller Container vector<Type> der Standardbibliothek int main(){ // Hauptfunktion int startTime = time(NULL); // Starten der Zeitmessung const int Anz_Pendel = 25; // Anzahl der zu berechnenden Pendel-Simulationen double v_phi_0 = 8.2; // Anfangsgeschwindigkeitsbereich (kleinster Wert) der Pendelsimulationen double beta = 0.0; // Anfangsgeschwindigkeitsbereich (kleinster Wert) der Pendelsimulationen double a = 0.0; // Untergrenze des Zeit-Intervalls double b = 50.0; // Obergrenze des Zeit-Intervalls int N_RK = 20000; // Anzahl der Gitter-Zeitpunkte des Runge-Kutta Ordnung vier Verfahrens unsigned int N_sim = 2000; // Anzahl der Gitter-Zeitpunkte der Simulation (Anzahl der ausgegebenen Punkte) double dt = (b - a)/N_sim; // Abstand dt zwischen den aequidistanten Punkten des Sim-Intervalls (h=dt) FILE *ausgabe; // Deklaration eines Files fuer die Ausgabedatei ausgabe = fopen("Pendel_Container.dat", "w+"); // Ergebnisse werden in die Datei "Pendel_Container.dat" geschrieben vector<Pendel> Kiste_Pendel; // Deklaration eines vector-Containers fuer die einzelnen Pendel-Objekte for (unsigned int n = 0; n < Anz_Pendel; ++n){ // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Pendel' if ((n % 5) == 0) { v_phi_0 = v_phi_0 * ( 1 + double(n)/55 ); beta = 0.0; } // Jeweils 5 Simulationen werden mit gleicher Anfangsgeschwindigkeit gemacht Kiste_Pendel.push_back( Pendel {a, 0.0, v_phi_0, 0.5, beta, N_sim} ); // Initialisierung der Pendel-Objekte beta = beta + 0.2; // Der Parameter beta der Reibung wird veraendert } // Beschreibung der in die Ausgabedatei "Pendel.dat" ausgegebenen Groessen fprintf(ausgabe, "%10s %20s %20s %20s %20s %20s %20s %20s \n", "# Index i", "t-Wert", "x(t), Pendel 1", "y(t), Pendel 1", "v_x(t), Pendel 1", "v_y(t), Pendel 1", "x(t), Pendel 2", "..."); // Zeitliche Simulation for(int i=0; i <= N_sim; ++i){ // for-Schleife ueber die Zeitgitterpunkte der Simulation fprintf(ausgabe, "%10d %20.12f ",i, Kiste_Pendel[0].t); // Ausgabe des Zeit-Indexes und der aktuellen Zeit in die Ausgabedatei for (auto& n : Kiste_Pendel){ // Bereichsbasierte for-Schleife ueber die einzelnen Pendel-Objekte fprintf(ausgabe, "%20.12f %20.12f %20.12f %20.12f ",n.r[0], n.r[1], n.v[0], n.v[1]); // Ausgabe der Orts- und Geschwindigkeitswerte des Pendel-Objektes in die Ausgabedatei n.Gehe_Zeitschritt(dt, N_RK); // Aufruf der Funktion Gehe_Zeitschritt(...) } // Ende for-Schleife (Pendel-Objekte) fprintf(ausgabe, "\n"); // Neue Zeile in der Ausgabedatei } // Ende for-Schleife (Zeitgitterpunkte) fclose(ausgabe); // Schliessen der Ausgabedatei cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; // Terminalausgabe der benoetigten Zeit } // Ende main()-Funktion
Da die einzelnen 25 Pendelsimulationen voneinander unabhängig sind, können wir die Parallelisierung an der folgenden for-Schleife anwenden: for (auto& n : Kiste_Pendel){ ... }. Um es OpenMP einfacher zu machen, schreiben wir diese bereichsbasierte for-Schleife zunächst als eine Index-basierte for-Schleife um (for (int n=0; n < Kiste_Pendel.size(); ++n){ ... }). Das folgende C++ stellt eine parallele Version des Programms dar:
/* Sequenzielle Version * Beispiel einer, von der Basisklasse 'Ding' abgeleitete Klasse 'Pendel' * Zwei zusaetzliche private Instanzvariablen (l und beta) und die winkelspezifischen Anfangswerte (phi und v_phi) * wurden der abgeleiteten Pendel-Klasse hinzugefuegt * Der Konstruktor der Klasse Pendel erzeugt ein Objekt Ding und initialisiert seine eigenen Daten-Member * Von der Klasse Pendel wurde eine weitere Subklasse 'Pendel_math' abgeleitet, welche die lineare Näherung der Pendel-DGL benutzt (gedaempfter harmonischer Ozillator) * Direkte Ausgabe der berechneten Werte von 25 Pendel-Simulationen in eine Datei (Visualisierung mittels Pendel_Container.py) * Zusatzliche Zeitmessung zum späteren Vergleich mit der parallelen Variante */ #include "Pendel.hpp" // Pendel und Ding Klassen #include <iostream> // Ein- und Ausgabebibliothek #include <vector> // Sequenzieller Container vector<Type> der Standardbibliothek #include <omp.h> // OpenMP zum parallelen Rechnen int main(){ // Hauptfunktion int startTime = time(NULL); // Starten der Zeitmessung const int Anz_Pendel = 25; // Anzahl der zu berechnenden Pendel-Simulationen double v_phi_0 = 8.2; // Anfangsgeschwindigkeitsbereich (kleinster Wert) der Pendelsimulationen double beta = 0.0; // Anfangsgeschwindigkeitsbereich (kleinster Wert) der Pendelsimulationen double a = 0.0; // Untergrenze des Zeit-Intervalls double b = 50.0; // Obergrenze des Zeit-Intervalls int N_RK = 20000; // Anzahl der Gitter-Zeitpunkte des Runge-Kutta Ordnung vier Verfahrens unsigned int N_sim = 2000; // Anzahl der Gitter-Zeitpunkte der Simulation (Anzahl der ausgegebenen Punkte) double dt = (b - a)/N_sim; // Abstand dt zwischen den aequidistanten Punkten des Sim-Intervalls (h=dt) FILE *ausgabe; // Deklaration eines Files fuer die Ausgabedatei ausgabe = fopen("Pendel_Container.dat", "w+"); // Ergebnisse werden in die Datei "Pendel_Container.dat" geschrieben vector<Pendel> Kiste_Pendel; // Deklaration eines vector-Containers fuer die einzelnen Pendel-Objekte for (unsigned int n = 0; n < Anz_Pendel; ++n){ // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Pendel' if ((n % 5) == 0) { v_phi_0 = v_phi_0 * ( 1 + double(n)/55 ); beta = 0.0; } // Jeweils 5 Simulationen werden mit gleicher Anfangsgeschwindigkeit gemacht Kiste_Pendel.push_back( Pendel {a, 0.0, v_phi_0, 0.5, beta, N_sim} ); // Initialisierung der Pendel-Objekte beta = beta + 0.2; // Der Parameter beta der Reibung wird veraendert } // Beschreibung der in die Ausgabedatei "Pendel.dat" ausgegebenen Groessen fprintf(ausgabe, "%10s %20s %20s %20s %20s %20s %20s %20s \n", "# Index i", "t-Wert", "x(t), Pendel 1", "y(t), Pendel 1", "v_x(t), Pendel 1", "v_y(t), Pendel 1", "x(t), Pendel 2", "..."); // Zeitliche Simulation for(int i=0; i <= N_sim; ++i){ // for-Schleife ueber die Zeitgitterpunkte der Simulation fprintf(ausgabe, "%10d %20.12f ",i, Kiste_Pendel[0].t); // Ausgabe des Zeit-Indexes und der aktuellen Zeit in die Ausgabedatei for (int n=0; n < Kiste_Pendel.size(); ++n){ // for-Schleife ueber die einzelnen Pendel-Objekte fprintf(ausgabe, "%20.12f %20.12f %20.12f %20.12f ",Kiste_Pendel[n].r[0], Kiste_Pendel[n].r[1], Kiste_Pendel[n].v[0], Kiste_Pendel[n].v[1]); } // Ende for-Schleife (Pendel-Objekte) fprintf(ausgabe, "\n"); // Neue Zeile in der Ausgabedatei #pragma omp parallel for // Parallelisierungsanweisung fuer die folgende for-Schleife for (int n=0; n < Kiste_Pendel.size(); ++n){ // for-Schleife ueber die einzelnen Pendel-Objekte Kiste_Pendel[n].Gehe_Zeitschritt(dt, N_RK); // Aufruf der Funktion Gehe_Zeitschritt(...) } // Ende for-Schleife (Pendel-Objekte) } // Ende for-Schleife (Zeitgitterpunkte) fclose(ausgabe); // Schliessen der Ausgabedatei cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; // Terminalausgabe der benoetigten Zeit } // Ende main()-Funktion
Da die einzelnen Threads nun immer einem separaten Pendelobjekt zugeordnet sind, ist eine zusätzliche private Kennzeichnung von Variablen nicht nötig und man kann einfach mittels "#pragma omp parallel for" parallelisieren. Die nebenstehende Abbildung zeigt die Terminalausgabe des sequentiellen und parallelen Programms und die unteren Bilder zeigen den zeitlichen Verlauf der Auslastung der Prozessorkerne. Man erkennt eine deutliche Reduzierung der benötigten Zeit.
