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} \]

Geom_Reihe_sequentiell.cpp
/* 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
    long 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(long 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:

Geom_Reihe_parallel.cpp
/* 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
    long 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(long 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.

Geom_Reihe_parallel_1.cpp
/* 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
    long 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(long 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:

Pendel_Container_sequentiell.cpp
/* Sequentielle Version
 * 25 Pendel-Simulationen in einem <vector>-Container
 * Das physikalische und mathematische Pendel mit Reibung (Stokesscher Ansatz)
 * Berechnung der Lösung einer Differentialgleichung zweiter Ordnung
 * bzw. eines Systems von zwei Differentialgleichungen erster Ordnung
 * mittels der Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) Pendel
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Konstruktor: Pendel(Masse des Pendels, Länge des Pendels l, Reibung beta, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha)
 * bzw. Konstruktor: Pendel_math(...)
 * Direkte Ausgabe der berechneten Werte von 25 Pendel-Simulationen in eine Datei (Visualisierung mittels Pendel_Container.py)
 */
#include "Pendel.hpp" // Header-Datei des physikalischen und mathematischen Pendels mit Reibung (Stokesscher Ansatz)
#include <iostream>   // Standard Input- und Output Bibliothek in C, z.B. printf(...)
#include <vector>     // Sequenzieller Container vector<Type> der Standardbibliothek
using namespace std;

// Hauptprogramm
int main() {
    int startTime = time(NULL);       // Starten der Zeitmessung
    size_t Anz_Pendel = 25;           // Anzahl der zu berechnenden Pendel-Simulationen
    double m = 1.0;                   // Masse des Pendels
    double l = 0.5;                   // Länge des Pendels
    double beta = 0;                  // Reibungskoeffizient
    double a = 0;                     // Untergrenze des Zeit-Intervalls
    double b = 20;                    // Obergrenze des Intervalls
    int N = 1000000;                  // Anzahl der Punkte
    vector<double> alpha = {0.0,8.2}; // Anfangswerte

    FILE *ausgabe;                                 // Deklaration eines Files fuer die Ausgabedatei
    ausgabe = fopen("Pendel_Container.dat", "w+"); // Ergebnisse werden in die Datei "Pendel.dat" geschrieben

    vector<Pendel> Kiste_Pendel(Anz_Pendel);       // Deklaration eines vector-Containers fuer die einzelnen Pendel-Objekte

    for (size_t n = 0; n < Anz_Pendel; ++n){                                // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Pendel'
        if ((n % 5) == 0) { alpha[1] *= (1 + double(n)/55); beta = 0.0; }   // Jeweils 5 Simulationen werden mit gleicher Anfangsgeschwindigkeit gemacht
        Kiste_Pendel[n] = Pendel(m,l,beta,a,b,N,alpha);                     // Initialisierung der Pendel-Objekte
        Kiste_Pendel[n].solve();                                            // Methode des Runge-Kutta-Verfahrens ausführen
        beta += 0.2;                                                        // Der Parameter beta der Reibung wird veraendert
    }

    // Beschreibung der in die Ausgabedatei "Pendel_Container.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", "...");

    // Ausgabe der Werte
    for(size_t i=0; i  < Kiste_Pendel[0].t.size(); ++i){                 // for-Schleife ueber die Zeitgitterpunkte der Simulation
        if ((i % 500) == 0) {                                            // nur jeden 100-ten Wert ausgeben lassen
            fprintf(ausgabe, "%10ld %20.12f ",i,  Kiste_Pendel[0].t[i]); // 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(i)[0],n.r(i)[1],n.v(i)[0],n.v(i)[1]); // Ausgabe der Orts- und Geschwindigkeitswerte des Pendel-Objektes
            }                                                           // Ende for-Schleife (Pendel-Objekte)
            fprintf(ausgabe, "\n");                                     // Neue Zeile in der Ausgabedatei
        }                                                               // Ende if
    }                                                                   // Ende for-Schleife (Zeitgitterpunkte)
    fclose(ausgabe);                                                    // Schliessen der Ausgabedatei

    cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; // Terminalausgabe der benoetigten Zeit
}                                                                                        // Ende Hauptprogramm

Da die einzelnen 25 Pendelsimulationen voneinander unabhängig sind, können wir das Programm in einfacher Weise parallelisieren. Wir lagern dazu die zeitintensive Berechnung der Lösung der Differentialgleichung 'Kiste_Pendel[n].solve();' in eine separate for-Schleife aus, die wir dann parallel ausführen. Das folgende C++ stellt eine parallele Version des Programms dar:

Pendel_Container_parallel.cpp
/* parallele Version
 * 25 Pendel-Simulationen in einem <vector>-Container
 * Das physikalische und mathematische Pendel mit Reibung (Stokesscher Ansatz)
 * Berechnung der Lösung einer Differentialgleichung zweiter Ordnung
 * bzw. eines Systems von zwei Differentialgleichungen erster Ordnung
 * mittels der Runge-Kutta Ordnung vier Methode
 * Verfahren zur Loesung der DGL ist in einer abstrakten Basisklasse ausgelagert
 * Berechnung der Lösung innerhalb einer Methode (Klassenfunktion) Runge-Kutta-Verfahren
 * Überschreibung der virtuellen Funktion der Bewegungsgleichungen in der abgeleiteten Klasse (Subklasse) Pendel
 * Zeitentwicklung der fuer unterschiedliche t-Werte in [a,b]
 * Konstruktor: Pendel(Masse des Pendels, Länge des Pendels l, Reibung beta, Anfangszeit a, Endzeit b, Anzahl der Punkte N, Anfangswerte alpha)
 * bzw. Konstruktor: Pendel_math(...)
 * Direkte Ausgabe der berechneten Werte von 25 Pendel-Simulationen in eine Datei (Visualisierung mittels Pendel_Container.py)
 */
#include "Pendel.hpp" // Header-Datei des physikalischen und mathematischen Pendels mit Reibung (Stokesscher Ansatz)
#include <iostream>   // Standard Input- und Output Bibliothek in C, z.B. printf(...)
#include <vector>     // Sequenzieller Container vector<Type> der Standardbibliothek
#include <omp.h>      // OpenMP zum parallelen Rechnen
using namespace std;

// Hauptprogramm
int main() {
    int startTime = time(NULL);       // Starten der Zeitmessung
    size_t Anz_Pendel = 25;           // Anzahl der zu berechnenden Pendel-Simulationen
    double m = 1.0;                   // Masse des Pendels
    double l = 0.5;                   // Länge des Pendels
    double beta = 0;                  // Reibungskoeffizient
    double a = 0;                     // Untergrenze des Zeit-Intervalls
    double b = 20;                    // Obergrenze des Intervalls
    int N = 1000000;                  // Anzahl der Punkte
    vector<double> alpha = {0.0,8.2}; // Anfangswerte

    FILE *ausgabe;                                 // Deklaration eines Files fuer die Ausgabedatei
    ausgabe = fopen("Pendel_Container.dat", "w+"); // Ergebnisse werden in die Datei "Pendel.dat" geschrieben

    vector<Pendel> Kiste_Pendel(Anz_Pendel);       // Deklaration eines vector-Containers fuer die einzelnen Pendel-Objekte

    for (size_t n = 0; n < Anz_Pendel; ++n){                                // for-Schleife zum Auffuellen des Containers mit Elementen vom Typ 'Pendel'
        if ((n % 5) == 0) { alpha[1] *= (1 + double(n)/55); beta = 0.0; }   // Jeweils 5 Simulationen werden mit gleicher Anfangsgeschwindigkeit gemacht
        Kiste_Pendel[n] = Pendel(m,l,beta,a,b,N,alpha);                     // Initialisierung der Pendel-Objekte
        beta += 0.2;                                                        // Der Parameter beta der Reibung wird veraendert
    }

    #pragma omp parallel for                                 // Parallelisierungsanweisung für die folgende, zeitintensive for-Schleife
    for (size_t n = 0; n < Anz_Pendel; ++n){                 // for-Schleife: Parallele Berechnung der Zeitentwicklung der einzelnen Pendel
        Kiste_Pendel[n].solve();                             // Methode des Runge-Kutta-Verfahrens ausführen
    }

    // Beschreibung der in die Ausgabedatei "Pendel_Container.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", "...");

    // Ausgabe der Werte
    for(size_t i=0; i  < Kiste_Pendel[0].t.size(); ++i){                 // for-Schleife ueber die Zeitgitterpunkte der Simulation
        if ((i % 500) == 0) {                                            // nur jeden 100-ten Wert ausgeben lassen
            fprintf(ausgabe, "%10ld %20.12f ",i,  Kiste_Pendel[0].t[i]); // 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(i)[0],n.r(i)[1],n.v(i)[0],n.v(i)[1]); // Ausgabe der Orts- und Geschwindigkeitswerte des Pendel-Objektes
            }                                                           // Ende for-Schleife (Pendel-Objekte)
            fprintf(ausgabe, "\n");                                     // Neue Zeile in der Ausgabedatei
        }                                                               // Ende if
    }                                                                   // Ende for-Schleife (Zeitgitterpunkte)
    fclose(ausgabe);                                                    // Schliessen der Ausgabedatei

    cout << "Das Programm benoetigte: " << time(NULL)-startTime << " Sekunden." << endl; // Terminalausgabe der benoetigten Zeit
}                                                                                        // Ende Hauptprogramm

Da die einzelnen Threads 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.