Dieses Semester (Sommersemester 2022) haben wir zum ersten mal im Rahmen der Veranstaltung Einführung in die Programmierung für Bio- und Chemieingenieure (kurz, EiP) eine Vorlesung zum Thema Effizienz von Code gehalten. Besonders teilenswert finde ich die Einfachheit, in der man heutzutage die Anzahl der Cores eines Rechnes ausnutzen kann, indem man For-Schleifen parallelisiert. Neben MATLAB, die Sprache, die in EiP gelehrt wird, bieten auch Sprachen wie C++ und Rust gute Möglichkeiten parallele For-Schleifen mit wenig Aufwand umzusetzen.

In diesem Beitrag betrachten wir ein Beispielproblem aus der Chemie: Das Berechnen der mittleren Temperaturen nach dem Zusammenmischen verschiedener Stoffe. Darauf aufbauend betrachten wir wie eine parfor-Schleife in MATLAB das Problem lösen kann. Umsetzungen in den systemnahen Programmiersprachen Rust und C++ werden wir uns gegeben falls in einem späteren Beitrag ansehen.

Beispielproblem

Die Mischtemperatur von mehreren Stoffen kann man mit der Richmannschen Mischungsregel berechnen, die von der folgenden Formel beschrieben wird:

$$T_m = \frac{\sum m_i \cdot c_i \cdot T_i}{\sum m_i \cdot c_i}$$

Als Eingabe benötigen wir die Wärmekapazitäten $~c_i~$, sowie die Massen $~m_i~$ als auch die Temperaturen $~T_i~$. Dazu werden wir zwei $~KxN~$ Matrizen anlegen. Hierbei entspricht N der Anzahl der Gemische, die berechnet werden und K die Anzahl der Stoffe die gemischt werden. Die erste Matrix T enthält die Temperaturen der K Gemisch-Komponenten für N Messpunkte und die zweite Matrix M enthält die Massen. Die Wärmekapazitäten speichern wir in einem K-elementigen Vektor. Im folgendem berechnen wir N mal die Temperatur $~T_m~$. Zuerst definieren wir unsere Daten Matrizen M und T, sowie die Wärmekapazitäten c mit dem folgenden Skript.

N = 1000;
%N = 1000000;
K = 4;
c = rand(1, K) * 0.5 + 2;
M = rand(N, K) * 30 + 1;
T = rand(N, K) * 50 + 5;

Wir führen also 1000 Berechnungen für Gemische mit vier Komponenten durch und legen dafür 4 zufällige Wärmekapazitäten im Bereich 2 bis 2.5 an, sowie 1000×4 Massen im Bereich 1 bis 31 und 1000×4 Temperaturen im Bereich 5-55 an.

Implementierung

Im folgendem verwenden wir Schleifen um mit den zufällig erzeugten Matrizen zu rechnen obwohl in MATLAB die Matrix-Operatoren sehr gut geeignet sind um das Problem ohne den Einsatz von Schleifen zu lösen. In produktiven Code wären Matrix-Operatoren die empfohlene Umsetzung. Eine mögliche Implementierung zur Berechnung der richmannsche Mischungstemperatur ist in dem folgendem Skript definiert:

results = zeros(N,1);
tic
for i=1:N 
%parfor i=1:N
    denominator = M(i,:) .* c(:);
    nominator = denominator(:) .* T(i, :);
    results (i) = sum(nominator) / sum(denominator);
end
fprintf("Duration with n=%i: %g seconds.\n", N, toc);

In Zeile 1 legen wir Speicher zum ablegen der Ergebnisse an. Danach starten wir mittels tic eine Zeitmessung. In Zeile 3/4 verwenden wir eine for- bzw. parfor-Schleife um die n Gemische zu berechnen. Im Schleifenrumpf wird die richmannsche Mischungsregel berechnet. Dazu wird zuerst der Nenner und dann der Zähler mittels elementweisen Matrixoperatoren berechnet um dann in Zeile 7 das Ergebnis zu berechnen. In Zeile 9 wird Ausgegeben wie lange die Berechnung gedauert hat.

Durch das Keyword parfor wird MATLAB veranlasst eine Datendekomposition durchzuführen. Intern werden die n Berechnungen in X Abschnitte unterteilt. Hierbei entspricht X häufig der Anzahl der Kerne in einem Rechner und jeder Abschnitt berechnet in etwa n/X Gemische. Dabei kann MATLAB manchmal vollautomatisch die Unterteilung vornehmen und kümmert sich um die nötige Kommunikation mit einem Threadpool, so dass der Programmierer damit keine großen Aufwände hat. Dazu müssen für eine parfor-Schleife bestimmte Bedingungen zutreffen: Es muss über Zählvariablen iteriert werden, Datenabhängigkeiten im Schleifen Rumpf, also z.B. ein Zugriff auf vorherigen Index i-1 sowie eine Verschachtelung von parfor-Schleifen ist nicht erlaubt. Weitereführende Informationen finden sich hier.

Zeitmessungen

Nun führen wir ein Zeitmessungen durch. Zuerst starten wir mit n=1000:

Duration with n=1000: 0.0061538 seconds. (for-Schleife)
Duration with n=1000: 0.105125 seconds. (parfor-Schleife)

Das ist unerwartet. Die parallele Variante benötigt mehr Zeit als die sequentielle. Erklärt werden kann dies mit der Initialisierungsarbeit, die MATLAB leisten muss um die Daten auf verschiedene Threads zu verschieben. Wir wiederholen das Experiment mit n=1000000:

Duration with n=1000000: 4.46868 seconds. (for-Schleife)
Duration with n=1000000: 1.16883 seconds. (parfor-Schleife)

Erst mit der größeren Datenmenge macht eine Parallelisierung des wenig aufwendigen Problems Sinn. Wenn die Operationen im Schleifen-Rumpf aufwändiger wären, als ein paar Additionen und Multiplikationen würde man den vorteilhaften Effekt der Parallelisierung schon bei einem kleinerem n beobachten können.

Zusammenfassung

Durch das Schlüsselwort parfor ist eine Schnelle Umsetzung von Parallelsierung auf Grundlage von Datendekomposition in MATLAB möglich. Wir haben gesehen, dass der Initialisierungsaufwand für kleine Probleme die Anwendung von parfor nicht rechtfertigen, aber für größere Probleme die Parallelisierung einen signfikanten Einfluss hat. Das schöne ist, dass der Einsatz von parallelen For-Schleifen oft ohne große Anpassungen möglich ist.