76 High Performance Computing
Tabelle 9.2:Installations- und Performanceparameter der Plattformen sgi Altix 3700 (Aster) und cray Opteron Cluster (Strider)
Aster Strider
cpu
Typ Intel Itanium 2 amd
Taktfrequenz (GHz) 1,3 2
Speicher (GB) 2 2
Performance (GFlops) 5,2 4
System
Anzahl Knoten 5 128
Anzahlcpus 416 256
Speicher (GB) 832 512
Performance (TFlops) 2,2 1,024
Verbindung IXS 16 GB s−1 Myrinet 2000
Zusammenfassend wird im Rahmen dieser Arbeit den Anforderungen der L¨osungsmethoden entsprechend folgende Parallelisierung der Analysealgorithmen vorgenommen:
• Umsetzung der direkten L¨osungsmethode auf einem ccnumaSystem
Dazu wird ein Knoten des sgi Altix 3700 Systems (Aster) als Rechenplattform mit gemeinsamem Cache genutzt.
• Umsetzung der iterativen L¨osungsmethode auf einemdmpSystem
Dazu wird dercrayOpteron Cluster (Strider) f¨ur die Parallelisierung verwendet.
Parallelisierung der direkten L¨osungsmethode 77
aj
aj aj
Aufteilung der Beobachtungen n
in bBlöcke Start Daten einlesen
Berechnung vonN undb
Ende CPU 1
CPU 1 CPU 2CPU 2
...
CPU CPUb N xˆ= −1
pmax
aj
aj aj
Aufteilung der Beobachtungen n
in bBlöcke Start Daten einlesen
Berechnung vonN undb
Ende CPU 1
CPU 1 CPU 2CPU 2
...
CPU CPUb N xˆ= −1
pmax
Abbildung 9.4: Parallelisierung der direkten L¨osungsmethode auf einemsmp oder ccnumaSystem
diejenige zum Aufbau der Designmatrix ¨ubersteigen). Mit diesem Hintergrund ist die durchgef¨uhrte Paral-lelisierung der direkten L¨osungsmethode in Abb. 9.4 graphisch dargestellt.
Der Analysealgorithmus wird bis zum Aufbau der Designmatrix sequentiell durchgef¨uhrt. Danach teilen sich pmax cpus den Aufwand zur Berechnung des Normalgleichungssystems. Die Inversion des Systems sowie die Ausgabe der Sch¨atzung ˆxerfolgt wiederum sequentiell. Der parallele Schritt ist in Abb. 9.5 genauer heraus gearbeitet. Die einzelnen Datenbl¨ocke der Dimension nb werden nach und nach abgearbeitet. Innerhalb eines Blockesi erfolgt die Aufteilung der Beobachtungen auf die beteiligten cpus, das heißt die Parallelisierung greift bei der Verarbeitung der Beobachtungen innerhalb eines jeden Blocks. Dabei berechnet jedecpuseine individuellen jp Zeilen der Designmatrix, es gilt Ppmax
p=1 jp = nb. Die einzelnen cpus greifen solange auf die Beobachtungen innerhalb des Blockesi zu, bis s¨amtliche Zeilen vonAi berechnet sind. Danach erfolgt die Aktualisierung der globalen FelderN=N+ATi Aiundb=b+ATiyi. Diese werden durch die Verwendung numerischer Bibliotheken basierend auf dem lapack/blas Standard (Anderson et al. 1999) automatisch parallel ausgef¨uhrt. Die parallele Umgebung endet mit der kompletten Aufstellung des Normalgleichungssy-stems (es ist im ¨ubrigen m¨oglich deren Inversion in die parallele Umgebung mit einzuschließen).
Die Implementierung gem¨aß Abb. 9.5 ist jedoch mit einem großen Nachteil verbunden. Die einzelnencpus greifen ohne fest definierte Reihenfolge auf die Beobachtungen zu, da deren Schnelligkeit und Auslastung durchaus variiert. Entsprechend werden diese in zuf¨alliger Reihenfolge prozessiert und demnach spiegelt der Aufbau der Designmatrix nicht zwingend den zeitlichen Verlauf der Beobachtungen wieder. Diese Tatsache wird dann zur Schwierigkeit, wenn eine fest vorgegebene Reihenfolge der Datenprozessierung vorausgesetzt wird. Dies ist z.B. der Fall beim Auftreten von zeitlichen Korrelationen zwischen den Beobachtungen. Diese werden im Allgemeinen in Form der Dekorrelation, vgl. Kapitel 6.1, ber¨ucksichtigt. Die zeitliche Dekorrelati-on verlangt eine Prozessierung der Daten in zeitlich fortlaufender Reihenfolge. Außerdem muss gew¨ahrleistet sein, dass die bereits dekorrelierten Werte mit einer bestimmten Korrelationsl¨ange vorliegen, um die Dekor-relation der aktuellen Beobachtung voran zu treiben. Besch¨aftigen sich jedoch mehrere cpus zeitgleich in besagter Korrelationsl¨ange, so ist eine derartige Prozessierung nicht m¨oglich.
Letztlich ist es also erforderlich, die a priori gew¨unschte Reihenfolge der Prozessierung zu erzwingen. Dies kann prinzipiell auf zwei Arten erfolgen. Zum einen durch entsprechende Restriktionen, zum anderen durch eine anders zu w¨ahlende Parallelisierung. Die erste M¨oglichkeit kann durch Bedingungsabfragen realisiert werden und wirkt sich verheerend auf die Performance der parallelen Region aus — verst¨andlich, da die Prozessoren immer wieder gebremst werden und viel Leerlauf haben. Folglich f¨uhrt nur Variante zwei zum Erfolg. Dabei erfolgt die Verteilung dernBeobachtungen auf diepmaxProzessoren derart, dass eine a priori festgelegte Anzahl (bzw. Reihenfolge) von Beobachtungen jeder einzelnencpuzugeteilt wird. Die pn
max
Beob-78 High Performance Computing
aj
aj aj
Blocki=1
Zeile aj=A( ,1:u) Zeile aj=A( ,1: u) Zeile aj=A( ,1: u) private jp=1
Σjp=n/b jp=jp+1
N= N+ Ni = N+AiTAi ; b= b+ bi = b+AiTyi
i=b i=i+1
CPU 1
CPU 1 CPU 2CPU 2
...
CPU CPUnein
ja
nein
ja
• • •
pmax
Abbildung 9.5:Parallelisierungsschritt gem¨aß Abb. 9.4 — Methode 1
achtungen pro Prozessor werden autonom abgearbeitet. Die Parallelisierung setzt so ein, dass nicht innerhalb der Bl¨ocke parallelisiert wird, sondern diese selbst auf diecpus aufgeteilt werden. Abbildung 9.6 zeigt dieses Vorgehen. Jedecputeilt ihrepn
max Beobachtungen in Prozessierungsbl¨ocke ein. Damit entstehen einzig durch die Aufteilung der Daten auf diecpus Bereiche, welche f¨ur die Dekorrelation nicht ber¨ucksichtigt sind. Dies kann durch die ¨Uberlappung der Beobachtungen um mindestens die Korrelationsl¨ange (in der praktischen Anwendung h¨angt die Dimension der ¨Uberlappung von der Warmlaufphase des Dekorrelationsfilters ab) be-hoben werden.
aj
aj aj
Block i=1
Zeile aj=A( ,1:u) Zeile aj=A( ,1: u) Zeile aj=A( ,1: u)
j=n/b j=j+1
i= b/pmax
i=i+1
CPU 1
CPU 1 CPU 2CPU 2
...
CPU CPUnein
ja
nein
ja
• • •
Block i=b/pmax+1
Block i=(pmax-1)b/pmax+1
j=1 j=1 j=1
j=n/b j=n/b
i= 2b/pmax i= b
j=j+1 nein
ja
j=j+1 nein
ja
ja ja
i=i+1 nein
i=i+1 nein
N= N+ Ni = N+AiTAi ; b= b+ bi = b+AiTyi
pmax
Abbildung 9.6:Parallelisierungsschritt gem¨aß Abb. 9.4 — Methode 2
Das alternative Parallelisierungsvorgehen bringt jedoch eine Reihe von ¨Anderungen in der Implementierung mit sich. Jede cpu erh¨alt ihre lokale Designmatrix Ai(nb ×u), was speichertechnisch unproblematisch ist.
Jedoch gilt dies nicht f¨ur die Normalgleichungsmatrix. Diese soll weiterhin global definiert sein (und davon keine Kopien existieren).
Parallelisierung der direkten L¨osungsmethode 79
Somit stellt sich die Aufgabe, die lokale Designmatrix pro Prozessor in den globalen Variablen N und b gem¨aß ATA und ATy zu aktualisieren. Zwei M¨oglichkeiten wurden f¨ur dieses Vorgehen untersucht. (i ) Jedecpuschreibt nacheinander in die globalen Variablen. Dieses Vorgehen kommt einer Entparallelisierung der algebraischen Operationen gleich und ist damit enorm langsam. (ii) Die algebraischen Operationen werden derart parallel ausgef¨uhrt, dass jeder Prozessor auf einen limitierten Speicherbereich der globalen Variablen schreibt (balancierte Aufteilung). Durch Permutation des Zugriffbereichs schreibt jedecpuletztlich in die gesamten globalen Variablen. Damit es zu keiner ¨Uberlappung in der Berechnung kommt, m¨ussen vor jeder Permutation die Prozessoren synchronisiert werden. Durch dieses Vorgehen ist ein Parallelisierungsgrad
¨ahnlich dem Prinzip nach Abb. 9.5 erreicht. Allein die Synchronisation reduziert eine optimale Auslastung dercpus.
Numerische Ergebnisse. Tabelle 9.3 fasst die erzielten Ergebnisse der direkten L¨osungsmethode zusam-men. Die Berechnungen erfolgen auf einemsgiAltix 3700 System, wobei bis zu 16 Prozessoren parallel zum Einsatz kommen. Deutlich erkennbar nimmt die Berechnung vonN=ATAundb=Ay, das heißt die Auf-stellung des Normalgleichungssystems, den Großteil des rechentechnischen Aufwands ein. Der prozentuale Anteil an der Gesamtlaufzeit betr¨agt meist ¨uber 90%.
Im Vergleich dazu nimmt der Aufbau der Designmatrix A nur einen unwesentlichen Posten ein. Dies gilt teilweise auch f¨ur die Inversion des Normalgleichungssystems. Auffallend ist jedoch, dass die Inversionspro-zedur mit steigender Anzahl ancpus nicht wesentlich beschleunigt. Mit 16 Prozessoren ist gar eine deutliche Verlangsamung zu beobachten, was am zunehmenden Kommunikationsaufwand zwischen den einzelnen Pro-zessierungselementen begr¨undet liegt.
Der Speed-up bzw. die Effizienz der parallelen Implementierung nimmt f¨urL = 200 ordentliche Werte an.
F¨ur L= 100 gilt dies ebenfalls bis einschließlich achtcpus. Mit der Verwendung von 16 cpus sinktEp auf 50 % . Tats¨achlich wurde dieses Verhalten schon in Austen et al. (2006) beobachtet.
Tabelle 9.3:Laufzeitergebnisse der direkten L¨osungsmethode auf der sgi Altix 3700 (Aster),n= 518 400 Beobachtungen
L= 100 L= 200
Anzahl Prozessoren Anzahl Prozessoren
1 4 8 16 1 4 8 16
Wall time (h) 3,84 1,10 0,66 0,48 64,46 17,06 9,06 5,57
Dateninput (min) 0,52 0,48 0,50 1,03 0,68 0,50 0,53 0,95 Referenzfeldabzug (min) 4,03 1.00 0,53 0,33 150,05 37,77 19,02 9,60 AufbauA(min) 5,84 2,08 1,32 0,77 31,68 8,10 4,32 2,77 Anteil an wall time (%) 2,5 3,2 3,3 2,7 0,8 0,8 0,8 0,8 BerechnungN,b(h) 3,64 1,07 0,61 0,42 61,43 15,71 8,05 4,46
Anteil an wall time (%) 94,8 97,3 92,4 87,5 93,8 92,1 88,8 80,0 Inversion (min) 1,52 0,63 0,67 1,35 59,38 34,67 37,08 53,35
Anteil an wall time (%) 0,7 1,0 1,7 4,7 1,5 3,4 6,8 15,6
Speed-upSp 1 3,49 5,82 8,00 1 3,84 7,22 11,76
EffizienzEp(%) 100 87 73 50 100 96 90 73
80 High Performance Computing