Effiziente parallele Implementierung eines expliziten Euler-Verfahrens für Grafikprozessoren durch Diamant-Tiling

Effiziente parallele Implementierung eines expliziten Euler-Verfahrens für Grafikprozessoren durch Diamant-Tiling Julien Kulbe Bayreuth Reports on Parallel and Distributed Systems No. 3, March 2012 University of Bayreuth Department of Mathematics, Physics and Computer Science Applied Computer Science 2 – Parallel and Distributed Systems 95440 Bayreuth Germany Phone: +49 921 55 7701 Fax: +49 921 55 7702 E-Mail: brp[email protected]yreuth.de Eziente parallele Implementierung eines expliziten Euler-Verfahrens für Grakprozessoren durch Diamant-Tiling Masterarbeit Universität Bayreuth Lehrstuhl für Angwandte Informatik 2 Betreuer: Prof. Dr. Th. Rauber, Dr. M. Korch und Dr. C. Scholtes Julien Kulbe 20. März 2012 Ich erkläre hiermit, • dass ich die vorliegende Arbeit ohne fremde Hilfe und ohne Verwendung anderer als der angegebenen Hilfsmittel verfasst habe, • dass ich sämtliche verwendeten Quellen erwähnt und gemäÿ gängigen wissenschaftlichen Zitierregeln korrekt zitiert habe. Inhaltsverzeichnis 1 Einführung 1 1.1 Motivation................................ 1 1.2 Ziele................................... 1 1.3 AufbauderArbeit ........................... 2 2 OpenCL 3 2.1 Aufbau eines OpenCL Programms . . . . . . . . . . . . . . . . . . 4 2.2 Optimierung eines OpenCL Programms . . . . . . . . . . . . . . . 6 2.2.1 LoopUnrolling......................... 6 2.2.2 Vektorisierung ......................... 7 2.2.3 Anzahl der Workitems und Workgroups . . . . . . . . . . . 7 2.2.4 Leseund Schreibzugrie auf den globalen Speicher . . . . . 8 2.2.5 Lokaler Speicher . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2.6 Barriers............................. 9 3 ODE-Systeme und Lösungsverfahren 11 3.1 Denition von ODE-Systemen . . . . . . . . . . . . . . . . . . . . . 11 3.2 Explizites Euler-Verfahren . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Ausgewählte Problemstellungen . . . . . . . . . . . . . . . . . . . . 14 4 CPU Implementierung des expliziten Euler-Verfahrens 17 4.1 Sequentielle Version . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.2 Pthreads parallelisierte Version . . . . . . . . . . . . . . . . . . . . 17 5 Lineare OpenCL Versionen 21 5.1 Unoptimierte lineare OpenCL-Version . . . . . . . . . . . . . . . . 21 5.2 Verbesserung des Speicherzugrismusters . . . . . . . . . . . . . . . 25 5.3 Vektorisierung der Problemfunktion . . . . . . . . . . . . . . . . . . 29 5.4 Probleme und Erkenntnisse bezüglich des linearen OpenCL-Verfahrens 30 6 Diamant-Tiling 31 6.1 Motivation................................ 31 6.2 Synchronisation............................. 31 6.3 Speicherverfahren............................ 32 6.4 Bestimmung der Diamantgröÿe . . . . . . . . . . . . . . . . . . . . 33 i Inhaltsverzeichnis 7 Zeilenweises Diamant-Tiling 35 7.1 Implementierung des Hostcodes . . . . . . . . . . . . . . . . . . . . 36 7.1.1 Vergröÿerung der Anzahl der Diamanten . . . . . . . . . . . 37 7.2 Implementierung des Kernels . . . . . . . . . . . . . . . . . . . . . 42 7.3 Laufzeit und Speedup des zeilenweisen Diamant-Tilings . . . . . . 45 7.4 Probleme und Erkenntnisse des Diamant-Tilings . . . . . . . . . . 51 8 Spaltenweises Diamant-Tiling 53 8.1 Implementierung des Hostcodes . . . . . . . . . . . . . . . . . . . . 53 8.2 Implementierung des Kernels . . . . . . . . . . . . . . . . . . . . . 54 8.3 Laufzeit des spaltenweisen Diamant-Tilings . . . . . . . . . . . . . 58 8.4 Probleme und Erkenntnisse des spaltenbasierten Diamant-Tilings . 58 9 Waben-Tiling 61 9.1 Motivation................................ 61 9.2 AufbauderWaben........................... 61 9.3 Ermitteln der optimalen Wabengröÿe . . . . . . . . . . . . . . . . . 63 9.4 Laufzeit und Speedup des Waben-Tilings . . . . . . . . . . . . . . . 66 9.5 Probleme und Erkenntnisse des Waben-Tilings . . . . . . . . . . . 66 10 Evaluierung 69 10.1 Evaluierung zum Einsatz von Grakprozessoren . . . . . . . . . . . 69 10.2 Evaluierung zum Optimieren von Algorithmen unter Grakprozessoren................................... 70 11 Zusammenfassung 73 ii 1 Einführung 1.1 Motivation Viele naturwissenschaftliche Problemstellungen lassen sich auf mathematische Beschreibungen durch Dierentialgleichungssysteme zurückführen. Das Lösen dieser Aufgaben ist einer der Schwerpunkte der numerischen Mathematik. Da diese Systeme sehr groÿ sein können und das Lösen sehr lange Zeit in Anspruch nimmt, ist das Parallelisieren und Optimieren von Algorithmen zum Lösen von Dierentialgleichungssystemen ein interessantes Thema. Die heutigen Grakprozessoren bieten durch eine Vielzahl an Multiprozessoren einen wesentlich höheren Grad an Parallelität als CPUs. Aus diesem Grund stellt gerade das Parallelisieren auf modernen Grakkartenarchitekturen eine Herausforderung dar. Grakprozessoren bieten unterschiedliche Granularitätsstufen zum Parallelisieren. Sie bestehen meist aus mehreren Multiprozessoren und diese wiederum aus mehreren Kernen. Jede Granularitätsstufe besitzt dabei Eigenschaften und Besonderheiten, die beachtet und analysiert werden müssen, um ihr Potential vollständig auszuschöpfen. 1.2 Ziele Das Ziel der Arbeit ist es, einen ezienten Algorithmus zum Lösen von Dierentialgleichungen für Grakkarten zu entwickeln. Dazu wird das Euler-Verfahren verwendet, da es von der Struktur sehr einfach ist und viele Möglichkeiten zum Optimieren bietet. Aber auch andere Verfahren lassen sich für die GPU portieren. Dabei muss immer abgewägt werden, welche Optimierungen für das Verfahren angewendet werden können. Zuerst soll eine Implementierung entwickelt werden, die sich am Algorithmus des Euler-Verfahrens für CPU orientiert. Da sich in dieser Version der Algorithmus am bestehenden Code orientiert, werden hier nur einige wichtige Optimierungen für Grakprozessoren angewendet. Aufwendige Optimierungen, wie zum Beispiel unterschiedliche Speicherhierarchien werden in dieser Version nicht betrachtet, da dafür groÿe Änderungen am bestehenden Code notwendig sind. Diese erste Version, lineare OpenCL-Version genannt, dient als Grundlage für alle weiteren Versionen. Das Ziel dieser Implementierung ist es, zu zeigen, dass bestehender CPU Code mit wenigen Änderungen ezient auf Grakkarten portiert werden kann. Dabei wird die Laufzeit der GPU-Version auf verschiedenen Hardware Architekturen ermittelt und mit der Laufzeit der bestehenden CPU-Version verglichen. Das nächste Ziel ist es, einen Algorithmus zu entwickeln, der möglichst viele Ei1 KAPITEL 2. OPENCL sein, wie die Anzahl der Recheneinheiten. Ist die Anzahl der Workgroups jedoch gröÿer als die der Recheneinheiten, muss auf eine gleichmäÿige Verteilung der Arbeit auf die unterschiedlichen Recheneinheiten geachtet werden. Es sollte also die Anzahl der Recheneinheiten der Grakkarte abgefragt werden, um dann die Anzahl der Workgroups auf die Anzahl der Recheneinheiten zu setzen, beziehungsweise auf ein ganzzahliges Vielfaches davon, soweit dies möglich ist. Die Anzahl der Workitems einer Workgroup sollte mindestens so groÿ sein, wie von einer Recheneinheit parallel bearbeiten werden kann. Auf AMD GPUs wird diese Gröÿe als wavefront und bei NVIDIA als warpsize bezeichnet. Auf heutigen Grakkarten beträgt sie 32 oder 64 Workitems. Jede Recheneinheit besitzt eine begrenzte Anzahl an Registern, die auf die einzelnen Workitems verteilt wird. Je gröÿer die Anzahl der Register pro Workitem ist, desto geringer ist die maximal mögliche Anzahl an Workitems pro Workgroup. Die maximale Anzahl der Woritems pro Workgroup kann durch clGetKernelWorkGroupInfo abgefragt werden. 2.2.4 Leseund Schreibzugrie auf den globalen Speicher Zum Austausch von Daten zwischen Host und Device kann der Host auf den globalen Speicher des Devices zugreifen. Der globale Speicher ist allerdings der langsamste Speicher des Devices und aus diesem Grund sollte versucht werden die Leseund Schreibzugrie auf den globalen Speicher zu minimieren. Beim Zugri auf den globalen Speicher sollte darauf geachtet werden, welches Zugrismuster am besten vom Device unterstützt wird. Die Grakkarte ist optimiert für 128-Bit Werte (= vier oat-Werte). Es ist also ratsam, jeweils 128 Bit am Stück zu lesen oder zu schreiben. Des Weiteren sollte angestrebt werden, dass alle Workitems gleichzeitig benachbarte Blöcke der Gröÿe 128-Bit laden. Das ideale Zugrismuster (nach [8, S. 583]) liegt vor, wenn jedes Workitem bei seinem Workitemindex beginnt oat4 (oder int4) Werte zu lesen. Dieser Anfangswert wird in den folgenden Iterationen um die Gesamtzahl der Workitems erhöht bis alle Daten aus dem Globalen Speicher gelesen wurden. Dies führt dazu, dass alle Workitems gleichzeitig einen groÿen Block von jeweils 128-bit Werten lesen (siehe Abbildung 2.2). uint id = globale Workitem Id uint ids = Anzahl globaler Workitems for ( idx = id ; idx < filterWidth ; idx += ids ) { 4 . . . } 2.2.5 Lokaler Speicher Um die Anzahl der Leseund Schreibzugrie des globalen Speichers zu optimieren, ist es sinnvoll, häug benutzte Daten zuerst in den lokalen Speicher zu laden und erst am Ende der Berechnung wieder in den globalen Speicher zu speichern. Jede Recheneinheit verfügt über einen lokalen Speicher. Dieser ist typischerweise 8 2.2. OPTIMIERUNG EINES OPENCL PROGRAMMS id=0 id=1 id=2 id=3 1. Iteration 2. Iteration id=0 id=1 id=2 id=3 ids 3. Iteration id=0 id=1 id=2 id=3 Abbildung 2.2: Zugrismuster des globalen Speichers wesentlich schneller als der globale Speicher. Es können jedoch nur die Workitems einer Workgroup auf den lokalen Speicher der Recheneinheit zugreifen, auf der sie abgearbeitet werden. Das stellt für Daten, die von mehreren Workgroups gebraucht werden, ein Problem dar. Um Daten zwischen verschiedenen lokalen Speichern auszutauschen, müssen sie über den globalen Speicher ausgetauscht werden. Dazu ist aber eine Synchronisierung der Workitems notwendig. 2.2.6 Barriers Barriers dienen zur Synchronisierung von Workitems. OpenCL bietet Barriers aber nur pro Workgroup an. Eine globale Barrier für alle Workitems ist durch die Verteilung der Workgroups auf verschiedene Recheneinheiten nicht praktikabel. Eine globale Barrier lässt sich realisieren, indem ein Kernel beendet wird und sich somit auch alle Workitems beenden. Der Host kann nun einen neuen Kernel starten und synchronisiert somit alle Workitems. Da dies relativ lange dauert, sollte versucht werden, die Anzahl der globalen Barriers zu minimieren und Sie durch lokale Barriers zu ersetzen. 9 3 ODE-Systeme und Lösungsverfahren 3.1 Denition von ODE-Systemen Viele naturwissenschaftliche Probleme lassen sich auf gewöhnliche Dierentialgleichungen ( ODE , ordinary dierential equation) mit gegebenem Anfangswert (IV, initial value) zurückführen. Solche Probleme werden als IVPs (initial value problems) bezeichnet. Ein IVP besteht zum einen aus Gleichungen, die das zeitliche Verhalten des Systems beschreiben (Gleichung 3.1) und zum anderen aus dem Anfangswert, der das System zum Zeitpunkt t0 beschreibt (Gleichung 3.2). Die Lösung des Problems besteht in der Bestimmung des Zustands des Systems zu einem bestimmten Zeitpunkt te> t0 (ausführlichere Informationen nden sich in [3] und [13]). Abbildung 3.1 zeigt ein exemplarisches ODE-Problem. Die durchgezogene Kurve zeigt den physikalisch exakten Verlauf, der jedoch bis auf den Anfangswert bei y0 unbekannt ist. Die gestrichelte Kurve zeigt den berechneten Verlauf eines Lösungsverfahrens für das gegebene ODE-Problem mit den berechneten Punkten yi zu den Zeitpunkten ti . y0(t) = f(t, y(t)) (3.1) y(t0) = y0 (3.2) Der Ausgangswert wird als n-stelliger Vektor modelliert. Aus diesem Ausgangswert werden iterativ mit einer gegebenen Schrittweite h neue Vektoren mit Werten berechnet, die den Zustand des System zum Zeitpunkt ti repräsentieren. Das explizite Euler-Verfahren ist eine Methode zur Lösung dieser Systeme. Dieses Verfahren ist eine einfaches Verfahren, um IVPs zu lösen. Es wird dafür eine fest vorgegebene Schrittweite h > 0 gewählt. Mit dieser Schrittweite werden dann die Zustände des Systems zu den Zeitpunkten ti=t0+i·h, i ∈N (3.3) berechnet. 11 KAPITEL 3. ODE-SYSTEME UND LÖSUNGSVERFAHREN y t t0t1t2te Physikalisch exakter Verlauf Numerisch berechneter Verlauf Abbildung 3.1: Berechnetes und exaktes Verhalten eines IVPs 3.2 Explizites Euler-Verfahren Listing 3.2 zeigt den Pseudocode des expliziten Euler-Verfahrens. H beschreibt dabei die Zeitdierenz zwischen t0 und te . In jedem Zeitschritt werden aus dem n-stelligen Vektor y _ cur für t die neuen Werte y _ new für den Zeitpunkt t+h berechnet. Die nächste Iteration verwendet dann y _ new als Ausgangswerte für den Zeitschritt t+h . In der Implementierung wird dies dadurch gelöst, dass y _ cur und y _ new vertauscht werden. Die Funktion ode_eval_comp beschreibt das spezi- sche Problem, welches gelöst werden soll. Sie berechnet für eine gegebene Komponente j des Vektors aus den aktuellen Werten die Steigung der Werte zum Zeitpunkt t. Diese Änderung, wird mit der Schrittweite h gewichtet, auf den aktuellen Wert addiert: yj,t+h=yj,t +h·ode _ eval _ comp(j, t, yt) (3.4) Daraus ergibt sich der neue Wert der Komponente j. Die Methode berechnet genau den Eintrag mit Index j der rechten Seite des ODE-Systems aus Gleichung 3.1. Diese neu berechneten Werte besitzen keine gegenseitigen Abhängigkeiten und können daher auch parallel zueinander berechnet werden. Die Berechnung der Lösung für ein spezisches Problem kann einen hohen Bedarf an Rechenaufwand darstellen, da zum einen die Anzahl der Komponenten sehr groÿ sein kann und zum anderen die Auswertung der Methode ode_eval_comp sehr komplex sein kann. Aufgrund des hohen Rechenaufwands und des hohen Grades an Parallelität bietet es sich an, ODEs auf GPUs zu berechnen. Abbildung 3.3 zeigt das prinzipielle Vorgehen dafür, wie aus den n Ausgangswerten iterativ die nächsten n Werte bis zum Endzeitpunkt te berechnet werden. Die schattierten Werte stellen dabei Einträge des Vektors dar, welche schon be12 3.2. EXPLIZITES EULER-VERFAHREN y_cur = Ausgangswerte_des_ODE_Systems ; for ( i = 0; i < H / h ; i++) { t = t_0 + i ∗ h; for ( j = 0; j < n; j++) { 5 y_new[ j ] = y_cur [ j ] + h ∗ ode_eval_comp( j , t , y_cur) ; swap(y_new, y_cur) ; } } print y_cur Abbildung 3.2: Pseudocode explizites Eulerverfahren rechnet wurden. Die weiÿen Einträge kennzeichnen dagegen noch nicht berechnete Werte. Zur Berechnung jedes Elements y _ new[j] , j je 1,...,n für den Zeitpunkt ti+1 werden bis zu 2·k+ 1 Ausgangswerte vom Zeitpunkt ti benötigt. Im allgemeinen Fall könnten alle Werte vom Zeitschritt ti gebraucht werden zur Berechnung von y _ new[j] . Bei vielen naturwissenschaftlichen Problemen lässt sich aber ein k bestimmen welches wesentlich kleiner als die Systemdimension n ist. Die Gröÿe k , die das Speicherzugrismuster und die Lokalität der Funktion f charakterisiert, ist die Zugrisdistanz. Die Zugrisdistanz gibt an, wie groÿ der Abstand vom zu berechnenden Wert mit Index j zum Index j+k des benötigten Wertes maximal ist. ti ti+1 te k n t0 t1 y_cur y_new k Abbildung 3.3: prinzipielles Vorgehen zur Berechnung 13 KAPITEL 3. ODE-SYSTEME UND LÖSUNGSVERFAHREN 3.3 Ausgewählte Problemstellungen Als Beispiele für simulierende Systeme wurden der Brüsselator (Bruss2d-Problem) und ein Modell für das Schwingen einer Saite (String-Problem) verwendet. Der Brüsselator ist ein Modell zur Beschreibung chemischer Oszillatoren. In Abbildung 3.4 ist das Brüsselator-Problem dargestellt. Das System besteht aus zwei verschiedenen Stoen mit den Konzentrationen A und B in einem zweidimensionalen Gitter. Das zweidimensionale Gitter besteht aus m∗m Gitterpunkten. Der n-stellige Vektor des resultierenden Systems repräsentiert die Konzentrationen der beiden Stoe in den m2 Gitterpunkten. Jedes Kästchen von Sto A und B steht dabei für eine Konzentration, die in dem n-stelligen Vektor abgespeichert wird. Um die Werte des zweidimensionalen Gitters des Brüsselators im Speicher abzubilden, werden die Werte zeilenweise jeweils abwechselnd für die Konzentration für Sto A und B im Speicher abgelegt. Zur Berechnung der neuen Konzentration von Sto A mit Punkt P = 22, werden die Werte der Konzentration von A von allen angrenzenden Gitterpunkten von P und von P selbst, sowie die Konzentration von Sto B im Punkt P benötigt. Beim Brüsselator mit einem m∗m Gitter dessen Werte zeilenweise abgelegt sind, entspricht die Zugrisdistanz einer vollständigen Zeile des Gitters. Sie beträgt 2m, da jeweils die Werte von Sto A und B pro Zeile abgelegt sind. Die Zugrisdistanz wächst demnach mit der Systemgröÿe. 4,4 5,3 5,4 5,5 6,4 m 20 56 56 51 57 63 52 64 58 Stoff B 20 21 5,3 28 23 24 m Stoff A 24 63 645,4 5,44,4 6,45251 5,3 16 5,5 2 * m 2 * m 4,2 21 2316 2858 57 59 22 Zugriffsdistanz 59 22 Abbildung 3.4: Abbildung des Brüsselators im Speicher Das String-Problem ([3, S. 27]) beschreibt das System einer schwingenden Saite, die an beiden Enden fest eingespannt ist. Die Saite wird durch äquidistante Abstände in Massepunkte unterteilt. Der n-stellige Vektor repräsentiert die Auslenkungen sowie die Geschwindigkeiten in den verschiedenen Massepunkten. Um die neue 14 3.3. AUSGEWÄHLTE PROBLEMSTELLUNGEN Auslenkung von Punkt P zu berechnen, wird die Geschwindigkeit von P selbst, der vorhergehenden Punkt P-1 und der nachfolgenden Punkt P+1 gebraucht. Die Zugrisdistanz beträgt dabei 3. Die Zugrisdistanz ist also hier unabhängig von der Gröÿe des Systems (siehe Abbildung 5.7). 00 1 2 3 4 5 678 9 10 5 5'4'43'3 3 Abbildung 3.5: String Problem 15 4 CPU Implementierung des expliziten Euler-Verfahrens 4.1 Sequentielle Version Als Ausgangspunkt der Optimierung steht eine Implementierung des expliziten Euler-Verfahrens von Dr. M. Korch [7] zur Verfügung. Listing 4.1 zeigt den zentralen Teil der Berechnung. Die Parameter der Funktion sind: die Startzeit zu der die Berechnung beginnen soll ( t0 ), die Anfangswerte ( y0 ) zum Startzeitpunkt, die Zeitspanne ( H =te−t0 ) über die berechnet werden soll, ein Zeiger auf das Feld, in dem der Lösungsvektor ( y ) eingetragen werden soll und die Schrittweite ( h ). Die Systemgröÿe n wird durch das zu lösende ODE-Problem vorgegeben ( ode_size ). Die Felder y0 und y besitzen jeweils n Einträge für die n Anfangswerte beziehungsweise Endwerte des Systems. Im Algorithmus wird zuerst ein Feld Y allokiert, in dem die Zwischenwerte der Berechnung gespeichert werden. Dieses Feld ist 2∗n groÿ, da es genügend Speicherplatz für die aktuellen Werte, sowie für die neu berechneten Werte des Systems bieten muss. Die Anzahl der Berechnungsschritte ( steps ) wird durch steps = Zeitspanne über die berechnet wird (H) Schrittweite (h) (4.1) errechnet. Danach werden die Ausgangswerte von y0 in das Feld Y kopiert, um die Berechnung zu initialisieren. In den einzelnen Berechnungsschritten wird jeweils das aktuelle ( Y_cur ) und das neue Feld ( Y_new ) alternierend bestimmt. Anschlieÿend werden für jeden einzelnen Wert des ODE-Systems die neuen Werte aus den aktuellen Werten gemäÿ der Formel 3.4 berechnet. 4.2 Pthreads parallelisierte Version Es werden in der Pthreads Version zwei Datenstrukturen verwendet, die über Argumente verfügbar sind. Einmal eine Datenstruktur die sich alle Threads teilen und eine für jeden Thread privat. In der gemeinsamen Datenstruktur stehen die Anzahl der Threads, die Startzeit, die Felder mit den Anfangs-, Endund Zwischenwerten und eine globale Pthreads-Barrier. Jeder Thread besitzt darüber hinaus noch eine private Datenstruktur für die Threadnummer und den Startbzw. Endindex der Komponenten, die er in jedem Zeitschritt berechnen soll. In der Berechnungsfunktion werden zuerst die Daten für die gemeinsamen und privaten Datenstrukturen 17 KAPITEL 5. LINEARE OPENCL VERSIONEN einen vergleichbaren Wert zu erhalten (Pthreads-Version mit 4 Threads, OpenCLVersion mit 512 lokalen Workitems und 30 Workgroups). Auf der x-Achse ist die Systemgröÿe aufgetragen und auf der y-Achse die entsprechende Laufzeit in Sekunden. Es wurden verschiedene Laufzeiten gemessen. Die Gesamtlaufzeit ist die komplette Laufzeit des Lösungsverfahrens mit Berechnung und Initialisierung. Die Kernellaufzeit ist dagegen die Laufzeit, die für die Berechnung gebraucht wurde. Es ist zu sehen, dass bei der Pthreads-Version die Gesamtund Kernellaufzeit nahezu identisch sind, da die Initialisierung bei Pthreads, die Erstellung der Threads, nicht lange dauert. Bei der OpenCL-Version ist dagegen die Initialisierungszeit wesentlich höher. Bei OpenCL muss zuerst die Grakkarte gefunden und initialisiert werden, der Kernel muss geladen und kompiliert werden und die Daten müssen zwischen Host und Device transferiert werden. Dieser Overhead zur Initialisierung ist für verschiedene Systemgröÿen nahezu konstant und deshalb nur bei sehr kleinen Systemen für die Gesamtlaufzeit signikant. Die Pthreads-Version ist für kleine Systeme schneller als die OpenCL-Version, da die Initialisierungszeit bei OpenCL hauptsächlich für die Laufzeit verantwortlich ist. Für groÿe Systeme, bei der die Initialisierungszeit keinen wesentlichen Beitrag zur Laufzeit mehr beisteuert, ist die OpenCL-Version doppelt so schnell wie die Pthreads-Version. Es wird deutlich, dass schon dieser einfache unoptimierte Kernel auf einer GPU zu einer Verbesserung der Laufzeit gegenüber einer CPU-Version führt. 0.01 0.1 1 10 100 10000 100000 1e+06 1e+07 Laufzeiten in Sekunden Systemgroesse n lin. OpenCL Gesamtlaufzeit lin. OpenCL Kernelzeit Pthreads Gesamtlaufzeit Pthreads Kernelzeit Abbildung 5.1: Laufzeit der linearen OpenCL-Version und der Pthreads-Version 24 5.2. VERBESSERUNG DES SPEICHERZUGRIFFSMUSTERS 5.2 Verbesserung des Speicherzugrismusters Zur Verbesserung des OpenCL-Kernels wurde zuerst das Zugrismuster auf den globalen Speicher verändert (siehe dazu Punkt 2.2.4). Bisher wurde für jedes Workitem beziehungsweise für jeden Thread ein Anfangsund Endindex berechnet. Dieses Verhalten ist jedoch für heutige CPU-Architekturen optimiert, die nicht einzelne Werte laden, sondern ganze Cachezeilen. Die meisten GPUs besitzen aber keine Caches und laden typischerweise 128-bit Blöcke. In Listing 5.3 ist der Kernel mit verbessertem Zugrismuster dargestellt. Jedes Workitem berechnet dabei immer vier oat-Werte / 128-bit Blöcke. Der Anfangsindex jedes Workitems ist dabei seine globale Workitem Id * 4. Somit werden in der ersten Iteration der Schleife die ersten 4∗Workitemanzahl Werte von y berechnet. Für die nächste Iteration wird der Anfangsindex um genau diese Anzahl von 4∗Workitemanzahl weitergezählt. Dies führt dann zu dem in Kapitel 2.2.4 beschriebenen optimalen Speicherzugrismuster von GPUs. In jeder Iteration werden dann vier aktuelle Werte gelesen und vier neue Werte geschrieben. Dies führt implizit zu einem Loop Unrolling (vgl. Kapitel 2.2.1). for ( int i = 4 ∗ id ; i <= ode size − 4; i+= 4 ∗ ids ) { 2 y_new[ i ] = y_cur [ i ] + h ∗ ode_eval_comp( i , t , y_cur) y_new[ i +1] = y_cur [ i +1] + h ∗ ode_eval_comp( i +1, t , y_cur) y_new[ i +2] = y_cur [ i +2] + h ∗ ode_eval_comp( i +2, t , y_cur) y_new[ i +3] = y_cur [ i +3] + h ∗ ode_eval_comp( i +3, t , y_cur) } Listing 5.3: optimiertes Speicherzugrismuster des linearen OpenCL-Verfahrens In Grak 5.2 ist die Kernellaufzeit der unoptimierten und der optimierten linearen OpenCL-Version dargestellt. Es wird deutlich, dass bei kleinen Systemgröÿen durch die Optimierung keine Verbesserung der Laufzeiten zu messen ist. Allerdings führt das Speicherzugrismuster schon bei mittleren Systemgröÿen zu einer erheblichen Verbesserung der Laufzeit. Bei der gröÿten gemessenen Systemgröÿe führt die Optimierung zu einem Speedup von 3,6 im Vergleich zur nicht optimierten Version. In Grak 5.3 sind für das String-Problem die Gesamtlaufzeiten der beiden linearen OpenCL-Versionen der Pthreads-Version gegenübergestellt. Es ist wieder zu sehen, dass die OpenCL-Versionen bei kleinen Systemgröÿen durch die Initialisierung schlechtere Laufzeiten haben als die CPU-Version. Bei gröÿeren Systemen sind jedoch die GPU-Versionen wesentlich schneller, die optimierte Version erreicht dabei eine Beschleunigung bis zum 7,4-fachen der CPU-Version. Um die Laufzeiten miteinander vergleichen zu können, ist es notwendig, die optimale Anzahl an lokalen und globalen Workitems für jeden Kernel zu bestim25 KAPITEL 5. LINEARE OPENCL VERSIONEN men. Die lokalen Workitems, also die Workitems pro Workgroup, werden durch die NVIDIA GTX 280 auf maximal 512 begrenzt. In Grak 5.4 ist der Speedup der Kernellaufzeiten gegenüber den lokalen Workitems aufgetragen. Es wurde dabei bei jeder Kurve eine feste Anzahl an Workgroups eingestellt und die lokalen Workitems von eins bis 512 variiert. Feste Anzahl pro Workgroup bedeutet, dass die Anzahl der globalen Workitems sich aus: globale _ Workitems =lokale _ Workitems ∗Workgroups (5.1) berechnet. Der Speedup jeder Kurve wurde dabei jeweils auf die Laufzeit mit einem globalen Workitem bezogen. Bei einer Workgroup ist zu Beginn ein linearer Anstieg des Speedups zu beobachten, der später dann immer weiter abacht. Das Maximum von 250 wird bei 512 Workitems, also der maximalen Anzahl, erreicht. Die optimale Anzahl an Workitems pro Workgroup ist hier die maximale Anzahl an Workitems pro Workgroup für die entsprechende GPU. Es stellt sich also die Frage warum dies so ist, wenn nur 8 Kerne pro Recheneinheit existieren. Eine Besonderheit dieses Kernels sind die vielen Speicherzugrie auf den globalen Speicher. Je mehr Workitems zur Verfügung stehen, desto besser kann die GPU die Latenzen der Speicherzugrie verstecken. Anders sieht es bei 30 Workgroups aus, also einer vollen Auslastung der GPU. Hier gibt es nur bis ca. 64 Workitems pro Workgroup einen starken Anstieg des Speedups bis auf das 1100-fache (entspricht etwa einem Speedup vom 35-fachen pro Workgroup). Danach bleibt der Speedup konstant. Dies liegt daran, dass bei 30 Workgroups viel mehr globale Workitems aktiv sind und auf den Speicher zugreifen. Bei etwa 30 Workgroups mit jeweils 64 Workitems, also 1920 globalen Workitems, ist die Bandbreite des globalen Speichers ausgelastet und es kommt zu keinem weiteren Speedup mehr. In Grak 5.5 wurden die Workitems pro Workgroup konstant gehalten und die Anzahl der Workgroups variiert. Bei allen Kurven ist anfangs wieder ein linearer Anstieg des Speedups zu sehen, der dann abacht und sich einer Obergrenze nähert. Wenn die Kurve abacht, ist wieder die maximale Bandbreite erreicht. Je gröÿer die Anzahl der lokalen Workitems ist, desto weniger Workgroups werden benötigt, um diesen maximalen Speedup zu erreichen. Für diesen Kernel sollte die Anzahl an lokalen und globalen Workitems möglichst hoch gewählt werden um den maximalen Speedup zu erhalten (zum Beispiel lokale Workitems: 512, Workgroups: 30 auf der GTX 280). 26 5.2. VERBESSERUNG DES SPEICHERZUGRIFFSMUSTERS 0.01 0.1 1 10 100 10000 100000 1e+06 1e+07 Laufzeit Systemgroesse n String - lin. OpenCL String - opt. Speicherzugrff Bruss2d - lin. OpenCL Bruss2d - opt. Speicherzugriff Abbildung 5.2: Laufzeiten der linearen OpenCL Implementierungen 0.1 1 10 100 10000 100000 1e+06 1e+07 Laufzeit Systemgroesse n lin. OpenCL lin. OpenCL (opt. Speicherzugriff) Pthreads Abbildung 5.3: Laufzeiten der Implementierungen 27 KAPITEL 5. LINEARE OPENCL VERSIONEN 0 200 400 600 800 1000 1200 0 100 200 300 400 500 600 Speedup lokale Workitems Workgroup=1 Workgroup=30 Abbildung 5.4: Speedup der Workgroups 100 200 300 400 500 600 700 800 900 1000 1100 0 10 20 30 40 50 60 Speedup Workgroup 128 lokale Workitems 256 lokale Workitems 512 lokale Workitems Abbildung 5.5: Speedup der lokalen Workitems 28 5.3. VEKTORISIERUNG DER PROBLEMFUNKTION 5.3 Vektorisierung der Problemfunktion Da sich die vier Berechnungszeilen in Kernel 5.3 sehr ähnlich sehen, liegt es nahe, diese zu vektorisieren (vgl. Kapitel 2.2.2). Die Schleife sieht dann so aus: for ( int i = id ; i < ODE_SIZE / 4; i += ids ) { y_new[ i ] = y_cur [ i ] + h ∗ ode_eval_comp4( i , t , y_cur) ; } Die normale Problemfunktion berechnet für einen gegebenen Index genau einen Rückgabewert. Wird die Schleife vektorisiert, muss auch die Problemfunktion angepasst werden, sodass für einen gegebenen Index vier Rückgabewerte berechnet werden. In Listing 5.6 und 5.7 sind die originale und die vektoriserte Problemfunktion für das String-Problem abgebildet. Daran ist ersichtlich, dass das Umschreiben nicht trivial ist, da alle Feldzugrie auf die vektorisierten Werte angepasst werden müssen. Bei komplexeren Funktionen sollte daher abgewägt werden, ob sich der Aufwand für das Umschreiben der Funktion im Vergleich zum Geschwindigkeitsgewinn lohnt. float ode_eval_comp (..) 2 { if ( i % 2 == 0) return y [ i +1]; return STRING_mod_K ∗ STRING_mod_K ∗ 7 (y [ i − 3] − 2 ∗ y [ i − 1] + y [ i +1]) ; } Abbildung 5.6: original String-Problem float4 ode_eval_comp4 ( . . ) 2 { float4 res ; res . x = y [ i ] . y ; res . z = y [ i ] .w; 7 res . y = STRING_mod_K ∗ STRING_mod_K ∗ (y [ i − 1].z − 2 ∗ y [ i ] . x + y [ i ] . z) ; res .w = STRING_mod_K ∗ STRING_mod_K ∗ (y [ i ] . x − 2 ∗ y [ i ] . z + y [ i +1].x) ; 12 return res ; } Abbildung 5.7: vektorisiertes String-Problem In Grak 5.8 wird die Kernellaufzeit der vektorisierten Version des String-Problems mit der bisherigen Version verglichen. Es ist zu sehen, dass das Vektorisie29 KAPITEL 5. LINEARE OPENCL VERSIONEN ren der Problemfunktion einen Speedup um das 2,5-fache zur nicht vektorisierten Funktion gebracht hat. Bei der Berechnung gröÿerer Systeme ist es sinnvoll, die Problemfunktion umzuschreiben. 0.01 0.1 1 10 100 10000 100000 1e+06 1e+07 Kernellaufzeit Systemgroesse n lin. OpenCL lin. OpenCL (opt. Speicherzugriff) lin. OpenCL (vek. Problemf.) Pthreads Abbildung 5.8: Laufzeit der vektorisierten Problemfunktion 5.4 Probleme und Erkenntnisse bezüglich des linearen OpenCL-Verfahrens Durch das Loop-Unrolling, das Verändern der Zugrisstruktur auf den globalen Speicher und das Anpassen der Problemfunktion konnten Laufzeitverbesserungen erzielt werden. Um aber weitere Optimierungsschritte anzugehen, bedarf es einer genaueren Analyse des Algorithmus, für weitere Laufzeitverbesserungen. Ein Problem ist, dass in jedem Kernelaufruf nur ein Zeitschritt berechnet wird. Dies führt dazu, dass das Synchronisieren des Hosts viel Zeit benötigt. Wie in Kapitel 2.2.6 beschrieben, sollte versucht werden ein Groÿteil dieser globalen Barriers durch lokale zu ersetzen. Weiterhin laufen bisher alle Speicheroperationen über den globalen Speicher, welcher der langsamste Speicher ist. Es sollte also versucht werden, einen Groÿteil der Berechnungsdaten vom globalen Speicher in den lokalen Speicher zu verlegen, um teure I/O-Speicherzugrie zu sparen. 30 6 Diamant-Tiling 6.1 Motivation Viele Anfangswertprobleme der ODE-Systeme besitzen eine beschränkte Zugrisdistanz. Das heiÿt, die Zugrisdistanz ist wesentlich kleiner als die Problemgröÿe n . Für die folgenden Algorithmen wird davon ausgegangen, dass das Problem eine beschränkte Zugrisdistanz aufweist, anderenfalls muss einer der vorherigen Algorithmen benutzt werden. Das Prinzip des Diamant-Tilings wurde in [12] benutzt, um Finite-DierenzenMethoden (FDTD, nite-dierence time-domain) auf Mehrkernprozessoren zu berechnen. Durch den Einsatz von Diamant-Tiling wurde die Speicherwiederverwendung eröht und es konnten dadurch Laufzeitverbesserungen erzielt werden. Es wurden dabei auch unterschiedliche Tileformen untersucht. Dabei stellte sich heraus, dass die Diamant-Tiles die beste Performance für Datanabhängigkeiten bei FDTD bieten. In der vorliegenden Arbeit wird darauf aufgebaut und das Diamant-Tiling in verschiedenen Implementierungen verglichen. In Kapitel 9 wird auch eine alternative Tileform untersucht. 6.2 Synchronisation Basierend auf der Zugrisdistanz acc _ dist kann der n-stellige Vektor in Blöcke gleicher Gröÿe mit jeweils block _ size ( ≥ Zugrisdistanz acc _ dist ) Werten unterteilt werden. Dies hat den Vorteil, dass die Werte in Block j zum Zeitschritt ti nur ihren eigenen Block j , den vorhergehenden Block j−1 und den nachfolgenden Block j+ 1 benötigen (siehe Abbildung 6.1). Damit kann das Problem behoben werden, dass pro Kernelaufruf nur ein Zeitschritt berechnet wird. Die Idee ist, jeder Recheneinheit eine bestimmte Anzahl (= dia_blocks ) an Blöcken zuzuteilen. Aus diesen dia _ blocks Blöcken zum Zeitpunkt ti können dia _ blocks −2 neue Blöcke für den Zeitschritt ti+1 berechnet werden. Es werden pro berechnetem Zeitschritt zwei Blöcke weniger, da die äuÿersten Blöcke nicht mit berechnet werden können, weil sie Abhängigkeiten zu Blöcken anderer Workgroups haben. Der Vorteil dabei ist, dass nach einem Zeitschritt die Workitems der Workgroup über eine lokale Barrier synchronisiert werden können. Erst wenn die Anzahl der Blöcke nur noch 2 beträgt und keine weiteren Blöcke mehr berechnet werden können, muss der Kernel beendet werden und es wird global synchronisiert. Aus diesem Muster ergibt sich die untere Hälfte eines Diamanten (vgl. Abbildung 6.1). 31 KAPITEL 6. DIAMANT-TILING Die Diamanten werden auf die Workgroups / Recheneinheiten verteilt. Dabei wird ein ähnliches Muster benutzt wie in Kapitel 2.2.4. Jede Workgroup beginnt mit dem Diamanten ihrer Workgroup Id (Anm: eindeutiger Index jeder Workgroup), für jeden weiteren Diamanten wird der Iterationszähler um die Anzahl der Workgroups weitergezählt, bis alle Diamanten der Iteration berechnet wurden (in Abbildung 6.2 wurden die Diamanten auf zwei Workgroups verteilt). Die Berechnungen einer Workgroup sind unabhängig von denen anderer Workgroups und werden parallel zu diesen ausgeführt. Das Berechnen von mehreren Diamanten einer Workgroup erfolgt sequentiell. Erst wenn alle Workgroups ihre Diamanten berechnet haben, kann die nächste Iteration gestartet werden. Nachdem alle unteren Hälften berechnet wurden, wird deshalb eine globale Barrier benötigt, um die Konsistenz der berechneten Daten sicherzustellen. dia_blocks block_size Abbildung 6.1: Aufbau eines Diamanten 6.3 Speicherverfahren Da die Berechnung jedes Diamanten nur auf einer Recheneinheit erfolgt, können die Daten für den Diamant im lokalen Speicher der entsprechenden Workgroup abgelegt werden. Aus den Abbildungen 6.1 und 6.2 ist ersichtlich, dass die nächsten Diamanten die zwei äuÿersten Blöcke (schraerte Blöcke) jedes Zeitschrittes der unteren Hälfte für ihre Berechnung der oberen Hälfte des Diamanten benötigen. In jedem Zeitschritt der unteren Hälfte müssen also die zwei äuÿersten Blöcke wieder in den globalen Speicher geschrieben werden. In Abbildung 6.3 ist zu sehen, dass es durch das Vertauschen von y_new und y_cur in jedem Zeitschritt möglich ist, alle Blöcke des Diamanten zu speichern, ohne noch benötigte Werte zu überschreiben. In der Abbildung werden alle geraden 32 6.4. BESTIMMUNG DER DIAMANTGRÖßE Zeitschritte (grob schraert) in dem vorderen Teil des 2n Werte groÿen Vektors y geschrieben und die ungeraden Zeitschritte (fein schraert) im hinteren Teil des Vektors. Bei der Berechnung der oberen Hälfte des Diamanten müssen die zwei äuÿersten Blöcke, die von den vorhergehenden Diamanten berechnet wurden, geladen werden. Die Daten werden dabei vom globalen Speicher in den lokalen Speicher kopiert. Die Berechnungen der Zeitschritte innerhalb des Diamanten erfolgt nur noch auf Daten im lokalen Speicher. Die Berechnungszeit innerhalb des Diamanten verkürzt sich, da die Zugriszeiten auf den lokalen Speicher geringer sind als auf den globalen Speicher und nur noch auf Daten im lokalen Speicher zugegrien wird. Es entsteht aber auch ein Mehraufwand, da die Daten erst zwischen den verschiedenen Speichern kopiert werden müssen. 1212 2 12 12 1 2 1 2 1 2 1 2 1 2 1 1 1 2 1 2 1 2 1 2 Abbildung 6.2: Vollständiges Diamant-Schema, Aufteilung der Diamanten auf (2 Workgroups) 6.4 Bestimmung der Diamantgröÿe Eine ideale Aufteilung wäre es, die n Werte in so viele Diamanten wie es Recheneinheiten gibt zu unterteilen, damit jede Workgroup einen Diamant berechnet. Damit ergibt sich für die Anzahl der Blöcke pro Diamant: dia _ blocks =n block _ dist ∗Recheneinheiten (6.1) Um Fallunterscheidungen und degenerierte Diamanten zu verhindern, wird die Anzahl der dia _ blocks Blöcke immer auf die nächstgröÿere gerade Zahl erhöht, da für die Diamanten, wie sie in Abbildung 6.1 dargestellt sind, immer eine gerade Anzahl an Blöcken benötigt wird. Ein anderes wichtiges Entscheidungskriterium für die Blockanzahl pro Diamant ist die Gröÿe des lokalen Speichers. Ist dieser 33 KAPITEL 7. ZEILENWEISES DIAMANT-TILING mult dia_blocks local_size 1 690 22.144 2 340 10.944 Als Resultat erhält man folgendes: in geraden Iterationen berechnen 30 Workgroups zwei Diamanten mit jeweils 340 Blöcken in ungeraden Iterationen berechnen 29 Workgroups zwei Diamanten, zwei Workgroups berechnen nur einen Diamanten Laufzeiten der drei Strategien In Abbildung 7.4 ist die Kernellaufzeit mit den unterschiedlichen Strategien dargestellt für verschiedene Systemgröÿen n. Bei einer Systemgröÿe von über 60.000 reicht der vorhandene lokale Speicherplatz nicht mehr für alle Diamanten und einige Workgroups müssen mehrere Diamanten berechnen. Den gröÿten Laufzeitsprung macht dabei die zweite Strategie, da hier nur eine Workgroup zwei Diamanten berechnen muss und alle anderen Workgroups warten müssen, bis der erste seine Arbeit beendet hat. Bei der ersten und dritten Strategie müssen alle Workgroups zwei Diamanten berechnen, die aber nur die halbe Gröÿe haben. Bei einer weiteren Erhöhung der Systemgröÿe bis etwa 120.000 kommt es zu keiner weiteren Vergröÿerung der Laufzeit der zweiten Strategie, da schon die erste Workgroup zwei Diamanten voller Gröÿe berechnet und damit die höchste Laufzeit aufweist. Bei den anderen beiden Varianten kommt es zu einem gleichmäÿigen Anstieg der Laufzeit. Es ist zu sehen, dass die Sprünge bei der Erhöhung der Diamantanzahl signikant für die Laufzeit sind. Die Vergröÿerung der Blockanzahl pro Diamant trägt jedoch weniger zur Laufzeit bei. Meistens wird durch die dritte Strategie die minimale Laufzeit für verschiedene Systemgröÿen erzielt. Mit den ermittelten Werten für die charakteristischen Gröÿen des Kernels kann der Kernel initialisiert werden. Dem Kernel wird zusätzlich zu den vorherigen Argumente noch die Blockgröÿe block _ size , die Zugrisdistanz und die Anzahl der Blöcke dia _ blocks übergeben. Danach folgt die Schleife über die Integrationsschritte, die die einzelnen Kernelaufrufe startet. Da jeder Diamant so viele Zeilen hoch hoch wie er Blöcke breit ist, berechnet jeder Diamant auch dia _ blocks Zeitschritte. Die nächste Iteration startet dann genau (dia _ blocks/2) ∗h später als der vorherige. (siehe Listing 7.5). 40 7.1. IMPLEMENTIERUNG DES HOSTCODES 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 0.12 0.13 0.14 50000 60000 70000 80000 90000 100000 110000 120000 130000 Kernellaufzeit Systemgroesse n 1. Strategie 2. Strategie 3. Strategie Abbildung 7.4: unterschiedliche Diamantverteilungsstrategien . . . int steps = H / h; float t = t_0 ; int max_iter = (2 ∗ steps ) / dia_blocks ; 5 setze Kernelargument ( block_size ) setze Kernelargument ( dia_blocks ) for ( iter = 0; iter <= max_iter ; ++iter ) { 10 setze Kernelargument ( ite r ) setze Kernelargument ( t ) starte Kernel ( global , lokal ) warte auf Beendigung des Kernels 15 t = iter ∗ dia_tiles / 2 ∗ h_f ; } Abbildung 7.5: Kernel starten im Hostcode 41 KAPITEL 7. ZEILENWEISES DIAMANT-TILING 7.2 Implementierung des Kernels Der Ablauf des Kernels ist in Listing 7.6 schematisch dargestellt. Der Vektor y enthält die Ausgangsdaten, die vom Host zum Device kopiert wurden. Der Vektor y _ l ist ein Speicherbereich im lokalen Speicher der Recheneinheit in dem die Workgroups die Berechnungen durchführen. Jede Workgroup wird auf einer Recheneinheit ausgeführt und besitzt deshalb einen lokalen Speicherbereich, der unabhängig von den anderen Workgroups ist. Das lokale Feld jeder Workgroup ist dabei anfangs nicht vorinitialisiert. def __kernel solve ( . . . , __global float4 ∗ y , 3 __local float4 ∗ y_l) {diamanten = Anzahl der Diamanten pro Zeile steps = dia_tiles / 2 − 1; for (group = 0; group < diamanten ; group+=workgroups) { 8 // obere Diamanthä l f t e for ( step = steps ; step >= 0; step −− ){ lade Werte aus y nach y_l berechne neue Werte in y_l } 13 // untere Diamanthä l f t e for ( step = 1; step <= steps ; step++) { berechne neue Werte in y_l speichere Werte von y_l nach y } 18 barrier } } Abbildung 7.6: Pseudocode des zeilenweisen Diamant-Verfahrens Die Anzahl der pro Zeitschritt zu berechnenden Diamanten wurde zuvor im Host bestimmt. Jede Workgroup besitzt eine Schleife über die Anzahl der Diamanten, die sie zu berechnen hat. Da für einen Diamanten der gesamte lokale Speicher benutzt wird, müssen die Diamanten auf einer Recheneinehit nacheinander berechnet werden. Die Diamanten auf verschiedenen Recheneinheiten werden aber parallel bearbeitet. Die Variable steps gibt die Höhe bzw. die Anzahl der Zeitschritte pro Diamanthälfte an. step = 0 entspricht dem Zeitschritt in der Mitte des Diamanten, bei dem die maximale Anzahl an Blöcken zu berechnen ist. Jede andere Stufe step des Diamanten besitzt 2∗(steps −step) Blöcke. Zunächst wird die obere Hälfte des Diamanten berechnet (siehe Abbildung 7.7 von step 5 bis 0). In dieser Phase werden in jedem Zeitschritt 4 Blöcke (schraf- erte gelbe Blöcke) vom globalen Speicher in den lokalen Speicher kopiert. Die 42 7.2. IMPLEMENTIERUNG DES KERNELS Berechnung der neuen Blöcke (schraerte blaue Blöcke ) erfolgt danach nur noch im lokalen Speicher. Der lokale Speicher ist wieder in zwei Teilbereiche für y _ new und y _ cur unterteilt, die nach jedem Zeitschritt vertauscht werden. In Abbildung 7.7 ist zu sehen, dass in step = 5 in das linke Feld (= y _ cur ) die Werte vom globalen Speicher geladen werden und auf der rechten Seite des Feldes ( y _ new ) werden die neu berechneten Werte abgelegt. Im nächsten Zeitschritt ( step = 4 ) wird dann das rechte Feld zu y _ cur und die nächsten Werte vom globalen Speicher werden neben die zuvor berechneten Werten geladen. Die neu berechneten Werte werden dann im linken Feld abgespeichert und überschreiben die zuvor geladenen Werte, da diese nicht mehr gebraucht werden. globaler Speicher lokaler Speicher Diamantschema step 0 1 2 3 4 5 2 1 3 dia_blocks =12 step = 5 step = 4 step = 3 step = 0 step = 1 step = 2 step = 5 4 5 Abbildung 7.7: Speicherverfahren des zeilenweisen Diamant-Tilings 43 KAPITEL 7. ZEILENWEISES DIAMANT-TILING Das Laden, Speichern und Berechnen von Werten erfolgt wieder nach dem Schema aus Kapitel 2.2.4. Es wird vor jedem Laden, Speichern und Berechnen der Startund Endindex der Werte berechnet und von allen Workitems parallel bearbeitet. id = lokale Workitem Id n = Anzahl der lokalen Workitems offset = Startadresse des Diamanten im globalen Speicher block_end = Endadresse des Diamanten im globalen Speicher 5 local_offset = Offset des Diamanten im lokalen Speicher . . . int begin = o ffset + step ∗ block_size; int end = block_end − step ∗ block_size; 10 for ( i = begin + id ; i < end ; i += n) { local_i = i − offset + local_offset ; y_new[ local_i ] . x = y_cur [ local_i ] . x + h ∗ ode_eval_comp (4 ∗ i , t , &y_cur [ local_offset − offset ]) ; . . . 15 } Abbildung 7.8: Berechnungen beim Diamant-Verfahren im lokalen Speicher Nachdem die obere Hälfte berechnet wurde, kann die untere Hälfte berechnet werden (in Abbildung 7.7 unten von step 1 bis 5). In dieser Phase der Berechnung eines Diamanten sind schon alle benötigten Daten im lokalen Speicher vorhanden und es können sofort die nächsten Zeitschritte berechnet werden ohne Daten neu aus dem globalen Speicher anzufordern. Nach dem Berechnen jeder Zeile müssen die äuÿersten vier Blöcke wieder in den globalen Speicher geschrieben werden, da sie zur Berechnung der nächsten Diamanten benötigt werden. Nach der Abarbeitung des Diamanten besteht der komplette globale Vektor aus neu berechneten Werten und diese können für den nächsten Diamanten wieder als Ausgangswerte benutzt werden. (Anm: In Abbildung 7.7 besteht der globale Bereich nicht nur aus neu berechneten Werten nach step = 5 . Dies liegt daran, dass es nur ein Ausschnitt des globalen Speichers ist und die noch gelben Blöcke von anderen Diamanten berechnet werden.) 44 7.3. LAUFZEIT UND SPEEDUP DES ZEILENWEISEN DIAMANT-TILINGS 7.3 Laufzeit und Speedup des zeilenweisen Diamant-Tilings Testsystem Pthreads/GTX280 GTX580 CPU: Prozessor 2x AMD Opteron 270 Intel E5530 Anzahl Kerne 2 4 Taktung 2.0 GHz 2.4 GHz GPU: Grakkarte NVIDIA GTX 280 NVIDIA GTX 580 Recheneinheiten 30 16 Anzahl Kerne 240 512 Kerne pro Recheneinheit 8 32 globaler Speicher 1023 MB 1536 MB lokaler Speicher 16 kB 48 kB L2 Cache / 768 kB Betriebssystem Linux Kernel 2.6.37 Compiler gcc 4.6.2 Im Folgenden werden die gemessenen Speedups der OpenCL-Version mit Diamant-Tiling betrachtet. Es wurde für alle Messungen das String-Problem mit einer Systemgröÿe von 80.000 verwendet. In Abbildung 7.9 wird der Speedup in Abhängigkeit von der Anzahl der Workgroups mit 128 und 256 lokalen Workitems pro Workgroup dargestellt. Hierbei sollte beachtet werden, dass wegen dem wesentlich komplexeren Kernelcode nun mehr Register als in der linearen OpenCL-Version benötigt werden. Da jede Workgroup nur eine bestimmte Anzahl an Registern für alle Workitems besitzt, ist die maximale Anzahl lokaler Workitems nur noch 256. Für eine höhere Anzahl an Workitems kann ein identisches Verhalten beobachtet werden. Es gibt einen linearen Speedup bis zu 30 Workgroups, da bei 30 Workgroups alle 30 Recheneinheiten komplett ausgelastet sind. Bei über 30 Workgroups gibt es einen starken Abfall des Speedups, da es jetzt eine ungleichmäÿige Verteilung der Diamanten/Workgroups auf die Recheneinheiten gibt. Bei 60 Workgroups gibt es wieder eine Lastbalancierung und es ergibt sich der gleiche Speedup wie bei 30 Workgroups. Die Anzahl der Workgroups sollte in der Diamant-Version also möglichst ein Vielfaches der Recheneinheiten sein. In Grak 7.10 ist der Speedup in Abhängigkeit von der Anzahl lokaler Workitems dargestellt. Es ist wieder anfangs ein linearer Anstieg des Speedups festzustellen, der sich ab 125 lokalen Workitems einer Konstanten nähert. Der maximale Speedup von 750 liegt bei etwa 160 lokalen Workitems, es werden also wieder wesentlich 45 KAPITEL 7. ZEILENWEISES DIAMANT-TILING mehr lokale Workitems als Kerne pro Recheneinheit benötigt, um den maximalen Speedup zu erreichen. Dies lässt darauf schlieÿen, dass wieder Speicherbandbreiten den Speedup begrenzen. Die Anzahl der Workitems sollte also im Idealfall auf über 100 lokale Workitems gesetzt werden. In Abbildung 7.11 sind die Laufzeiten der Diamant-Versionen im Vergleich zu den vorherigen Algorithmen aufgetragen. Es werden zusätzlich zu dem bereits vorgestellten Diamant-Verfahren noch zwei weitere betrachtet. Ein Verfahren benutzt lediglich globalen Speicher und keinen lokalen Speicher und bei einem anderen wurde zusätzlich die Problemfunktion vektorisiert. Vergleichbar sind daher das Diamant-Verfahren mit dem linearen OpenCL-Verfahren und die beiden Verfahren mit der vektorisierten Problemfunktion. Die Diamant-Version erreicht dabei einen Speedup von 1,75-fache, der linearen OpenCL Version. Das heiÿt, dass die Verwendung des lokalen Speichers und der Einsatz lokaler Barriers die Laufeit des Kernels fast halbiert. Bei den Versionen mit der vektorisierten Problemfunktion ist der Speedup durch das Diamant-Verfahren jedoch wesentlich geringer. Bei kleinen Systemgröÿen wird noch ein Speedup von 1,5 erreicht, bei gröÿeren Systemen nur von 1,05. Die Ursache dafür ist in Abbildung 7.12 zu erkennen. In Abbildung 7.12 ist die normierte Kernellaufzeit gegenüber der Systemgröÿe aufgetragen. Die normierte Kernellaufzeit berechnet sich aus: normierte Kernellaufzeit = Kernellaufzeit Systemgröÿe n (7.5) Sie stellt die durchschnittliche Berechnungszeit für ein einzelnes Element dar. Bei beiden Diamant-Versionen ist zu sehen, dass mit wachsender Systemgröÿe die normierte Laufzeit geringer wird (bis zu einer Systemgröÿe von 60.000). Dies liegt daran, dass bei gröÿeren Systemen auch die einzelnen Diamanten gröÿer werden und damit der lokale Speicher besser ausgenutzt werden kann und es weniger globale Barriers gibt. Ab einer Systemgröÿe von 60.000 müssen aber mehrere Diamanten von einer Workgroup berechnet werden (siehe Abbildung 7.4). Deshalb steigt die Laufzeit sprunghaft wieder an. Bei 60.000, also einer vollen Ausnutzung des lokalen Speichers, wird ein Minimum der normierten Laufzeit erzielt, da bei einer vollen Ausnutzung des lokalen Speichers der höchste Speedup des Diamant-Verfahrens erreicht wird. Jede weitere Vergröÿerung des Systems bringt keinen weiteren Speedup mehr für jedes einzelne Element. Dieses Speedupmaximum wird immer wieder bei Vielfachen von 60.000 erreicht. Der Speedup des Diamant-Verfahrens ist also primär abhängig von der Gröÿe des lokalen Speichers. Die lineare OpenCL-Version dagegen hat einen stetigen Speedup pro Einzelelement. Bei kleinen Systemen ist die normierte Laufzeit wesentlich schlechter als die der Diamant-Versionen, bei gröÿeren Systemen wird aber annähernd die Laufzeit der Diamant-Version erreicht. Der Speedup der linearen OpenCL-Version ist stark abhängig von der Anzahl der globalen Barriers. Da bei allen Systemgröÿen die Anzahl der globalen Barriers konstant ist, fällt die Ausführungszeit für die globalen Barriers für kleine Systeme wesentlich stärker ins Gewicht als bei gröÿeren Systemen, bei denen die Berechnung der Elemente signikant für die Laufzeit ist. 46 7.3. LAUFZEIT UND SPEEDUP DES ZEILENWEISEN DIAMANT-TILINGS 0 100 200 300 400 500 600 700 800 0 10 20 30 40 50 60 Speedup Workgroups 128 lokale Workitems 256 lokale Workitems Abbildung 7.9: Speedup der lokalen Workitems 1 10 100 1000 5 25 50 125 300 Speedup lokale Workitems 1 Workgroup 30 workgroups Abbildung 7.10: Speedup der Workgroups 47 KAPITEL 7. ZEILENWEISES DIAMANT-TILING 0.01 0.1 1 10 100 10000 100000 1e+06 1e+07 Kernellaufzeit Systemgroesse n lineare OpenCL Version lin. OpenCL Version (vek. Problemf.) Diamant Tiling (vek. Problemf.) Diamant Tiling (ohne lokalen Speicher) Diamant Tiling Abbildung 7.11: Laufzeit der Versionen 5e-07 1e-06 1.5e-06 2e-06 2.5e-06 3e-06 3.5e-06 4e-06 4.5e-06 5e-06 5.5e-06 6e-06 10000 100000 1e+06 1e+07 normierte Kernellaufzeit Systemgroesse n lin. OpenCL (vek. Problemf.) Diamant Tiling Diamant Tiling (vek. Problemf.) Abbildung 7.12: normierte Laufzeit der Versionen 48 7.3. LAUFZEIT UND SPEEDUP DES ZEILENWEISEN DIAMANT-TILINGS Wenn der Speedup der Diamant-Version primär von der Gröÿe des lokalen Speichers abhängt, dann müsste sich die Laufzeitdierenz der Diamant-Version zur linearen OpenCL-Version bei GPUs mit einem gröÿeren lokalen Speicher vergröÿern. Um dies zu untersuchen, wurden die Laufzeitmessungen auch auf einer NVIDIA GTX 580 durchgeführt, die 48kB lokalen Speicher besitzt, im Gegensatz zur GTX 280 mit nur 16kB. Dabei muss jedoch beachtet werden, dass die GTX 580 einen L2 Cache besitzt. Jede Speicheranfrage auf den globalen Speicher lädt daher nicht nur einen 32bit Wert, sondern 128Byte, also eine komplette Cachezeile. In Abbildung 7.13 ist die normierte Laufzeit bei verschiedenen Problemgröÿen zu sehen. Die Laufzeiten der Diamant-Versionen und der linearen Versionen sind nahezu identisch. Dies liegt daran, dass die lineare Version linear auf den globalen Speicher zugreift und daher stark vom Cache protiert. Die Diamant-Version lädt immer nur einzelne Blöcke vom globalen Speicher und kann daher weit weniger vom Cache protieren als die lineare Version. Da die Diamant-Version für GPUs ohne Cache optimiert wurde, wurde in Abbildung 7.14 der Cache durch deklarieren der Vektoren als volatile umgangen. volatile bewirkt, dass die GPU gezwungen wird, jeden Wert zu laden und damit nicht vom Cache protieren kann. Die Laufzeiten der Diamant-Versionen bleiben dabei nahezu konstant, da sie nicht vom Cache protieren. Die Laufzeiten der linearen Versionen steigen jedoch. Auÿerdem ist zu beobachten, dass die linearen Versionen vom Vektorisieren der Problemfunktion stark protieren und beim Diamant-Verfahren das Vektorisieren zu fast keinem Speedup führt. Dies lässt darauf schlieÿen, dass beim Diamant-Verfahren die Speicherbandbreite erschöpft ist und, dass daher die Vektorisierung der Problemfunktion nur noch zu einer geringen Beschleunigung führt. In Abbildung 7.15 wurden die Laufzeiten für das Brüsselator-Problem auf der GTX 580 gemessen. Das Brüsselator-Problem hat eine wesentlich gröÿere Zugrisdistanz als das String-Problem, welche auch mit der Gröÿe des Systems wächst. Dies hat zur Folge, dass die Blöcke der Diamanten wesentlich gröÿer sind. Für das Diamant-Tiling wird eine gewisse Mindestanzahl an Blöcken vorausgesetzt, damit es eektiv arbeiten kann. Durch die geringe lokale Speichergröÿe wurde daher auf der GTX 280 auf Messungen mit dem Brüsselator-Problem verzichtet. Die Abbildung 7.15 zeigt die normierten Laufzeiten für verschiedene Systemgröÿen auf der GTX 580. Ohne Cache ist die Diamant-Version wesentlich schneller als die lineare Version, mit Cache ist die lineare Version besser als die Diamant-Version, da sie sehr stark vom L2 Cache protiert. 49 KAPITEL 8. SPALTENWEISES DIAMANT-TILING def __kernel solve ( . . . , __global float4 ∗ y , __local float4 ∗ y_l) { 5 groups = Anzahl der Diamanten // = mult steps = dia_tiles ; for (group = 0; group < groups ; group+=workgroups ) { for ( diag = 0; diag < steps ; diag++){ for ( height = 0; height < steps /2; height++) { 10 if ( diag == 0 | | height == 0) lade Werte aus y in y_l berechne Werte in y_l if ( diag >= steps − 2 | | height == steps − 1) speichere Werte aus y_l in y 15 barrier } } } } Abbildung 8.3: Pseudocode des spaltenweisen Diamant-Verfahrens Betrachtet wird das String-Problem mit einer Zugrisdistanz von 4 (Anm: Zugrisdistanz von 3 wird auf das nächst gröÿere Velfache von 4 aufgerundet, durch die Vektorisierung mit oat4-Werten) die Berechnung wird von jeweils 512 lokalen Workitems durchgeführt daraus ergibt sich eine Blockgröÿe von: block_size = 512 Workitems * 4 Werte = 2048 Werten Die Anzahl der Werte die beim spaltenweisen Diamant-Tiling geladen/gespeichert werden müssen pro Zeile: unoptimiert: alle vier Blöcke werden vollständig kopiert → 4 Blöcke * 2048 Werte = 8192 Werte optimiert: zwei vollständige Blöcke und zwei Blöcke bis zur Zugrisdistanz → 2 Blöcke * 2048 Werte + 2 Blöcke * 4 Werte = 4104 Werte Im Vergleich dazu die Anzahl der Werte die beim zeilenweisen Diamant-Tiling geladen/gespeichert werden pro Zeile: → 4 Blöcke * 4 Werte = 16 Werte → Durch die Optimierung konnte der Bedarf an zu speichernden Werten fast halbiert werden. Der Bedarf beim zeilenweisen Diamant-Tiling ist jedoch insbesondere für kleine Zugrisdistanzen wesentlich geringer. 56 8.2. IMPLEMENTIERUNG DES KERNELS lokaler Speicher Diamantschema dia_blocks =6 3 B C D E 1 2 GF H JI LK654 13 121110987 A 1817 161514 B CA B C 1 B C1 E F E F3 B C1E F3JI B C1E F3JI7 B C1E F3JI7 D2 B C1E F3JI7 D24 B C1E F3JI7 D24I8 B K1E F3JI17 L24I18 JI15 I16 JI11 I12 diag = 0 diag = 1 diag = 5 Abbildung 8.4: Spaltenweises Diamant-Tiling 57 KAPITEL 8. SPALTENWEISES DIAMANT-TILING 8.3 Laufzeit des spaltenweisen Diamant-Tilings Die Laufzeitmessungen für das spaltenweise Verfahren wurden alle auf der GTX 580 durchgeführt. Da das Verfahren für groÿe Blöcke ausgelegt ist, ist zu erwarten, dass gerade das Brüsselator-Problem von der Optimierung protiert. Für das Brüsselator-Problem können aber nur auf der GTX 580 gröÿere Messungen durchgeführt werden, da bei der GTX 280 der kleinere lokale Speicher nur für geringe Systemgröÿen Messungen zulässt. In Abbildung 8.5 ist die normierte Laufzeit des String-Problems für die GTX 580 dargestellt. Es ist zu sehen, dass für kleine Systemgröÿen keine Messungen für das spaltenbasierte Verfahren existieren, da durch die Vergröÿerung der Blöcke nicht genügend Blöcke existieren und somit das Diamant-Verfahren nicht eektiv arbeiten kann. Ab einer Systemgröÿe von 100.000 wird der gesamte lokale Speicher benutzt und der maximale Speedup des Verfahrens kann erreicht werden. Die Laufzeit ist jedoch gröÿer als beim zeilenweisen Verfahren. Dies kann mit der Vergröÿerung der Blöcke begründet werden. Zum einen verringert dies die Zeitschritte pro Diamant und führt damit zu mehr globalen Barriers. Zum anderen erhöhen dies die Anzahl der Werte, die zwischen den Speichern kopiert werden müssen. In Abbildung 8.6 ist die normierte Laufzeit für das Brüsselator-Problem für die GTX 580 dargestellt. Durch die gröÿere Zugrisdistanz des Brüsselator-Problems müssen die Blöcke ab einer bestimmten Problemgröÿe nicht weiter vergröÿert werden (ab n= 125.000 ). Es ist zu sehen, dass das spaltenbasierte Verfahren auch hier kaum eine Laufzeitverbesserung erzielen kann. Das liegt vor allem daran, dass der Speedup durch die eektivere Speicherausnutzung und durch die höhere Anzahl von lokalen Barriers kompensiert wird. 8.4 Probleme und Erkenntnisse des spaltenbasierten Diamant-Tilings Es hat sich herausgestellt, dass trotz der eektiveren Speicherausnutzung des spaltenbasierten Diamant-Tilings keine Laufzeitverbesserungen zum zeilenweisen Diamant-Tiling erzielt werden konnten. Dies liegt daran, dass zum einen die gröÿeren Blöcke zu mehr Speicherund Ladeoperationen führen und zum anderen mehr lokale Barriers benötigt werden und damit die Workitems häuger unterbrochen werden. Das spaltenbasierte Verfahren bietet daher keine Laufzeitverbesserungen zum zeilenweisen Diamant-Verfahren. Das spaltenweise Verfahren bietet aber dahingehend einen Vorteil, dass sich dadurch Probleme mit groÿen Zugrisdistanzen lösen lassen, für die ansonsten der Speicher beim zeilenweisen Verfahren nicht ausreichen würde. 58 8.4. PROBLEME UND ERKENNTNISSE DES SPALTENBASIERTEN DIAMANT-TILINGS 0 5e-07 1e-06 1.5e-06 2e-06 2.5e-06 3e-06 3.5e-06 10000 100000 1e+06 1e+07 normierte Kernellaufzeit Systemgroesse n ode=STRING_MOD lin. OpenCL zeilenweises Diamantt. spaltenweises Diamantt. Abbildung 8.5: normierte Kernellaufzeit String 4e-07 6e-07 8e-07 1e-06 1.2e-06 1.4e-06 1.6e-06 1.8e-06 2e-06 0 50000 100000 150000 200000 250000 300000 350000 400000 450000 500000 normierte Kernellaufzeit Systemgroesse n ode=BRUSS2D-MIX lin. OpenCL zeilenweises Diamantt. spaltenweises Diamantt. Abbildung 8.6: normierte Kernellaufzeit Bruss2d 59 9 Waben-Tiling 9.1 Motivation Das Waben-Tiling ist eine Variation des zeilenweisen Diamant-Tilings. Durch die beiden bisher vorgestellten Versionen des Diamant-Tilings wurde deutlich, wie wichtig es ist, die Anzahl der Barriers zu minimieren und die Arbeit zwischen den Barriers zu maximieren. Beim spaltenweisen Verfahren gab es viele lokale Barriers und wenig Arbeit für jedes Workitem, was sich negativ auf die Laufzeit ausgewirkt hat. Beim zeilenweisen Verfahren war die Anzahl der Barriers gering, die Arbeit pro Zeile jedoch stark unterschiedlich. An der breitesten Stelle bei step = 0 (siehe Abbildung 7.7) ist die Arbeit maximal, an den oberen und unteren Rändern des Diamanten sind jedoch nur wenige Elemente pro Zeile zu berechnen. Gerade bei kleinen Zugrisdistanzen, bei denen die Blöcke nur aus wenigen Elementen bestehen, werden an den Rändern nicht alle Workitems beschäftigt. Das WabenTiling versucht genau dies zu verbessern. Beim Diamant-Tiling berechnete jeder Diamant mit dia _ blocks Blöcken genau dia _ blocks−1 Zeitschritte. Beim WabenTiling wird eine feste Anzahl an Schritten pro Diamant vorgegeben (= steps ) und damit werden die Enden jedes Diamanten abgeschnitten. Das Ergebnis ist, dass auch an den oberen und unteren Rändern noch genügend Arbeit für alle Workitems vorhanden ist. Ein anderer Faktor der dazu beiträgt, dass die lineare OpenCL Version sehr gute Laufzeiten erzielt, ist, dass der lineare Zugri auf groÿe zusammenhängende Bereiche des globalen Speichers, wie er in Kapitel 2.2.4 beschrieben wurde, für die Grakkarten optimal ist. Dieses Speicherzugrismuster kann aber nur in der linearen Version erreicht werden. In den Diamantversionen werden pro Zeile nur zwei eher kleine zusammenhängende Bereiche zwischen globalen und lokalen Speicher übertragen. Beim Waben-Tiling kann der Speicher an den Randzeilen wieder linear übertragen werden und das optimale Zugrismuster kann erreicht werden. 9.2 Aufbau der Waben Abbildung 9.1 zeigt das prinzipielle Schema beim Waben-Tiling. Das Abarbeiten und Berechnen der Diamanten durch die Workgroups erfolgt wie bei der zeilenweisen Diamant-Version. Jede Workgroup berechnet ihren Diamanten zeilenweise. Eine Besonderheit des Waben-Tilings ist, dass die Waben ineinander verschränkt sind. Da die verschiedenen Iterationen sich nun überlagern, kann Formel 6.1 nicht mehr zur Aufteilung der Waben auf die Workgroups benutzt werden. Die Gröÿe 61 KAPITEL 9. WABEN-TILING Abbildung 9.1: Waben-Tiling eines Blockes wird wie beim zeilenweisen Diamant-Verfahren auf das nächste Vielfache von vier der Zugrisdistanz gesetzt, damit die Blöcke so klein wie möglich werden, um unnötige Speicheroperationen zu vermeiden. Die Gesamtzahl aller Blöcke ( dia _ blocks _ sum ) berechnet sich aus der Problemgröÿe geteilt durch die Gröÿe eines Blockes: dia _ blocks _ sum =n/block _ size (9.1) 2 1 0 3 3 Wabe Wabe Zwischenwabe 0 1 2 2 steps = 3 dia_blocks dia_blocks - 2 * steps dia_blocks Abbildung 9.2: Wabe und Zwischenwaben mit steps = 3 Beim zeilenweisen Diamant-Tiling wurden alle Blöcke auf die Workgroups verteilt. Beim Waben Tiling müssen die Blöcke jedoch auf die Waben (gelbe Waben) und die Zwischenwaben (blaue Waben) aufgeteilt werden. Zur Bestimmung dieser Aufteilung betrachte man einen Zeitschritt, der step = 0 einer Wabe entspricht. In Abbildung refcomb2 entspricht dies step = 3 einer Zwischenwabe. Die Aufteilung 62 9.3. ERMITTELN DER OPTIMALEN WABENGRÖßE wird so vorgenommen, dass es in solch einem Zeitschritt bei workgroups Waben immer workgroups−1 Zwischenwaben gibt. Die maximale Anzahl an Blöcken pro Zeitschritt der Waben wird wieder als dia _ blocks bezeichnet. In einem Zeitschritt, in dem eine Wabe bei step 0 die maximale Anzahl an Blöcken pro Zeitschritt berechnet, haben die entsprechenden Zwischenwaben gerade die kleinste Anzahl an Blöcken (siehe Abbildung refcomb2). Die kleinste Anzahl der Blöcke ist abhängig von dem Zeitschritt nach dem eine Wabe abgeschnitten wird. In der Abbildung wird jede Wabe und Zwischenwabe nach dem Zeitschritt 3 (= steps) nach oben und unten abgeschnitten. Die Anzahl der Blöcke bei der maximalen Anzahl von Zeitschritten ist dia _ blocks −2∗steps . Die Gesamtzahl aller Blöcke muss auf die workgroups Waben und die workgroups −1 Zwischenwaben verteilt werden: dia _ blocks _ sum = Blöcke der Waben + Blöcke der Zwischenwaben dia _ blocks _ sum =workgroups ∗dia _ blocks + (workgroups −1) ∗(dia _ blocks −2∗steps) Durch Umformung lässt sich daraus die Anzahl der Blöcke pro Diamant für eine bestimmte Schrittweite ( steps ) berechnen: dia _ blocks _ sum = 2 ∗workgroups ∗dia _ blocks+ (workgroups −1) ∗(−2∗steps)−dia _ blocks dia _ blocks _ sum + 2 ∗steps∗(workgroups −1) = (2 ∗workgroups −1) ∗dia _ blocks dia _ blocks =dia _ blocks _ sum + 2 ∗steps ∗(workgroups −1) 2∗workgroups −1 9.3 Ermitteln der optimalen Wabengröÿe Wird die Gröÿe der Waben minimiert ( steps = 0 ), so wird ein ähnliches Verhalten wie beim linearen OpenCL-Verfahren erzielt. Es wird dabei das optimale Zugrismuster für alle Werte vom globalen Speicher erzielt, da alle Werte auf einmal geladen werden. Jedoch werden die Daten im lokalen Speicher nicht wiederverwendet. Werden die Waben vergröÿert, erhält man einen positiven Eekt, indem die Werte im lokalen Speicher nun wiederverwendet werden. Allerdings gibt es auch einen negativen Eekt, da der Zugri auf die globalen Daten nicht mehr immer dem optimalen Muster folgt. Werte am Rand müssen in jedem Schritt wieder nachgeladen werden, was zu einem nicht optimalen Zugrisverhalten auf den globalen Speicher führt. Wird die Wabengröÿe maximiert, ergibt sich daraus das Verhalten des zeilenweisen Diamant-Tilings. Durch Variation der Gröÿe steps ist es möglich 63 KAPITEL 9. WABEN-TILING einen optimalen Wert zu erreichen, der für die jeweilige Hardware den Kompromiss zwischen einer Wiederverwendung des lokalen Speichers und den optimierten Speicherzugrien auf den globalen Speicher ndet. Nachfolgend sind die Laufzeiten des String-Problems mit einer Systemgröÿe von 2.000.000 dargestellt. In Abbildung 9.3 ist die Laufzeit für die GTX 280 mit ca. 500 Blöcken pro Diamant in Abhängigkeit von der maximalen Schrittweite steps aufgetragen. Erwartungsgemäÿ wird bei kleinen Schrittweiten mit wachsender Schrittweite die Laufzeit geringer, da eine Wiederverwendung durch den lokalen Speicher erzielt wird. Ab einer Schrittweite von ca. 50 steigt die Laufzeit wieder an, da der unoptimierte Zugri auf den globalen Speicher signikant wird. Bei der minimalen Laufzeit ist das Verhältnis der Wiederverwendung des lokalen Speichers zum optimierten Speicherzugri auf den globalen Speicher optimal. Die Wiederverwendung wird durch die max. Schrittweite charakterisiert und die Menge der Daten, die zwischen den Speichern transferiert werden muss, wird durch die Anzahl der Blöcke charakterisiert. Beim Minimum der Laufzeit beträgt das Verhältnis der Blöcke pro Diamant zur Schrittweite: ∼500/40 =∼12.5 . In Grak 9.4 ist die Abhängigkeit von der maximalen Schrittweite auf der GTX 580 dargestellt. Durch ihren gröÿeren lokalen Speicher besteht hier jeder Diamant aus ca. 1500 Blöcken. Es ist zu sehen, dass das Minimum bei etwa 100 Schritten erreicht ist. Dies ist höher als bei der GTX 280, jedoch im Verhältnis zur Blockgröÿe ( 1500/100 =∼15 ) stellt sich ein Wert ähnlicher Gröÿenordnung ein. Das bei manchen Schrittweiten auftretende plötzliche Ansteigen der Laufzeiten ist dadurch zu erklären, dass durch Variation der maximalen Schrittweite auch die Anzahl an Blöcken und damit die Anzahl an Diamanten verändert wird und dies zu Laufzeitsprüngen führen kann. Als Beispiel des betrachtet man den Laufzeitsprung in Abbildung 9.4: bei einer maximalen Schrittweite von 100 auf 120 kommt es zu einem Laufzeitsprung: Die GTX 580 besitzt 16 Recheneinheiten In der folgenden Tabelle sind für die max. Schrittweiten die Anzahl der Blöcke pro Diamant angegeben, die Gesamtanzahl alle Diamanten und das Verhältnis wenn alle Diamanten auf 16 Recheneinheiten aufgeteilt werden: Schrittweite Blöcke Diamanten #Diamanten / #16 Rechene. 100 1526 176 11 110 1528 177 11.06 120 1530 178 11.1 Durch die Änderung der maximalen Schrittweite ändert sich auch die Anzahl der Blöcke pro Diamant und die Anzahl der Diamanten. Bei 100 ist die Anzahl gerade ein Vielfaches der Recheneinheiten, es werden also genau 11 Diamanten von jeder Workgroup berechnet. Ab Schrittweite 110 muss eine Recheneinheit 12 Diamanten berechnen, was einen Sprung in der Laufzeit verursacht. 64 9.3. ERMITTELN DER OPTIMALEN WABENGRÖßE 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 1 10 100 normierte Kernellaufzeit max. Schrittweite 128 lokale Workitems (vek. Problemfunktion) 128 lokale Workitems 192 lokale Workitems (vek. Problemfunktion) 192lokale Workitems Abbildung 9.3: Laufzeiten für verschiedene Maximalschrittweiten GTX 280 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 10 100 1000 normierte Kernellaufzeit max. Schrittweite 512 lok. Worki. (vek. Problemfunktion) 512 lokale Workitems 256 lok. Worki. (vek. Problemfunktion) 256 lokale Workitems Abbildung 9.4: Laufzeiten für verschiedene Maximalschrittweiten GTX 580 65 11 Zusammenfassung Das Ziel der Arbeit war es, das Euler-Verfahren auf Grakprozessoren zu portieren und dafür zu optimieren. Dies ist zum einem mit dem linearen OpenCLVerfahren und zum anderen mit den Diamant-/Waben-Tilings gelungen. Das lineare OpenCL-Verfahren ist zwar ein sehr einfaches Verfahren, welches nicht alle Optimierungen der GPU ausnutzt, aber es werden gute Laufzeiten erzielt, die wesentlich schneller sind als vergleichbare CPU-Versionen. Algorithmen die schon für die CPU parallelisiert wurden, zum Beispiel mit Pthreads oder OpenMP, lassen sich leicht für die GPU portieren. Auf der GPU können dann schon einige wichtige aber einfache Optimierungen zu sehr guten Laufzeiten führen, die sich nur schwer mit der CPU erreichen lassen. Das Diamant-/WabenTiling ist wesentlich komplexer, führt aber auch meist zu noch besseren Ergebnissen als das lineare OpenCL-Verfahren, da es mehr Vorteile der Grakprozessoren ausnutzen kann. Bei sehr zeitintensiven Berechnungen, bei denen es wichtig ist, die Laufzeit zu minimieren, wäre es also von Vorteil den bestehenden Code derart anzupassen, dass möglichst alle Eigenschaften und Besonderheiten der GPUs berücksichtigt werden. Dabei muss beachtet werden, dass sich einige Optimierungen gegensätzlich verhalten können. Wie beim Waben-Tiling gesehen, musste ein Kompromiss zwischen einem optimalen Zugrismuster und der Wiederverwendung des lokalen Speichers gefunden werden. Der Einsatz von Grakprozessoren zum Lösen von Dierentialgleichungen ist bei Problemen, die lange zum Berechnen benötigen, zu empfehlen. Das Initialisieren der GPU und das Übersetzen des Kernels benötigen allerdings Zeit, die sich über die Laufzeit der Berechnung amortisieren muss. OpenCL hat sich im Laufe der Arbeit als eine geeingete Programmierschnittstelle für Grakprozessoren herausgestellt, bei der jedoch anzumerken ist, dass sie speziell für Grakprozessoren entwickelt wurde. OpenCL-Kernel lassen sich zwar auch auf CPUs ausführen, jedoch sind viele Besonderheiten auf Prozessoren nicht vorhanden und führen dort nicht immer zu einem Laufzeitgewinn, wie zum Beispiel der Einsatz lokaler Barriers und lokaler Speicher. Optimierungen wie zum Beispiel das Speicherzugrismuster, verhalten sich auf CPU und GPU teilweise nachteilig, sodass spezielle Kernel für GPU und CPU entwickelt werden müssen, um den optimalen Speedup auf allen Devices zu erhalten. Warum ndet die GPU Programmierung trotz enormer Geschwindigkeitssteigerungen gegenüber Prozessoren noch keine breite Verwendung in kommerziellen Programmen? Ein Grund kann darin gesehen werden, dass bisher keine standardisierte Schnittstelle für alle Grakprozessoren zur Verfügung stand. Jeder Grakkartenhersteller bevorzugte seine eigene API (NVIDIA CUDA, ATI Stream). Mit 73 KAPITEL 11. ZUSAMMENFASSUNG OpenCL steht eine geeignete Schnittstelle für alle Grakkarten zur Verfügung, die den Einsatz von GPUs für allgemeine Berechnungen weiter verbreiten kann. Die jüngste Vergangenheit hat gezeigt, dass die Steigerung der Geschwindigkeit von Prozessoren zum Groÿteil durch die Erhöhung der Anzahl der Kerne und nicht mehr so sehr durch Erhöhung des Prozessortaktes erreicht wird. Deshalb ist zu erwarten, dass der Einsatz von Grakkarten zur Berechnung von zeitintensiven Prozessen in Zukunft zunimmt. Einen wichtigen Faktor werden dabei PC Spiele darstellen, da sie den Systemen durch immer neue Innovationen wie verbesserte Physik Modellierung und KI immer mehr Leistung abverlangen. Beispiele, bei dem der Geschwindigkeitsgewinn von GPUs für nicht Grak-bezogene Aufgaben benutzt wird, sind Spiele wie Batman: Arkham City und Batman: Arkham Asylum ([5],[4]), die die Physik des Spiels durch GPU-Physix ([11]) berechnen lassen. Aber auch Aufgaben des wissenschaftlichen Rechnens, nutzen immer mehr Grakprozessoren zum Berechnen. Daher beziehen schon drei der fünf schnellsten Supercomputer ([1]) einen Groÿteil ihrer Rechenkapazitäten aus Grakkarten. Sollten sich in der Industrie einheitliche Standards bei der GPU-Programmierung durchsetzen, mit APIs, die auch über mehrere Jahre hinweg beständig sind, so ist zu erwarten, dass die GPU-Programmierung in Zukunft immer breitere Verwendung ndet. 74 Literaturverzeichnis [1] TOP500 Supercomputer Site. Website: http: // www . top500 . org , besucht am 13.03.2012. [2] Amazon. Website: http: // www . amazon . de , besucht am 2.12.2011. [3] E. Hairer, S.P. Nørsett, and G. Wanner. Solving ordinary dierential equations: Nonsti problems . Springer series in computational mathematics. Springer, 1993. [4] PC Games Hardware. Website: http: // www . pcgameshardware . de / aid , 8503 63 / Batman-Arkham-City-Trailer-zeigt-GPU-Physx-Update-Bildvergle ich-mit-und-ohne-Effekte-in-hoher-Aufloesung / Action-Spiel / News / , besucht am 14.12.2011. [5] Eidos Interacrive. Website: urlhttp://www.batmanarkhamcity.org/, besucht am 14.12.2011. [6] Khronos OpenCL Working Group. The OpenCL Specication, version 1.1 , 2011. [7] Matthias Korch. Eziente Implementierung eingebetteter Runge-KuttaVerfahren durch Ausnutzung der Speicherzugrislokalität . Doctoral thesis, University of Bayreuth, December 2006. [8] Adam Lake. Game Programming Gems 8 . Course Technology, 2010. [9] Microsoft. DirectCompute. Website: http: // www . microsoftpdc . com / 2009 / P0 9-16 , besucht am 14.12.2011. [10] NVIDIA. Cuda. Website: http: // www . nvidia . de / object / cuda \ _home \ _new \ _de . html , besucht am 14.12.2011. [11] NVIDIA. GPU-Physix. Website: http: // www . geforce . com / Hardware / Tech nologies / physx , besucht am 14.12.2011. [12] Daniel A. Orozco and Guang R. Gao. Mapping the FDTD Application to Many-Core Chip Architectures. In Proceedings of the 2009 International Conference on Parallel Processing , ICPP '09, pages 309316, Washington, DC, USA, 2009. IEEE Computer Society. [13] K. Strehmel and R. Weiner. Numerik gewöhnlicher Dierentialgleichungen . Teubner Studienbücher. B.G. Teubner, 1995. 75 Literaturverzeichnis [14] Kane S. Yee. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media. IEEE Trans. Antennas and Propagation , pages 302307, 1966. 76 Zusammenfassung Die hier vorliegende Arbeit beschäftigt sich damit, das explizite Euler-Verfahren auf Grakprozessoren zu optimieren. Dabei werden die Speicherhierarchien, lokale Datenwiederverwendung, Ausnutzung der Speicherbandbreite der GPU und die Synchronisierung zwischen Host und Device genauer untersucht. Dabei werden zwei Implementierungen näher betrachtet, das Diamant-Tiling und das lineare Verfahren, da sie sich gut eignen um die Optimierungen genauer zu untersuchen. Es stellt sich dabei heraus, dass Optimierungen wie die lokale Datenwiederverwendung und der optimale Zugri auf den Speicher sich gegensätzlich verhalten. Ein Mischverfahren (das Waben-Tiling), dass dabei die Vorteile des linearen Verfahrens und des Diamant-Tilings vereint, führt daher zu den besten Laufzeiten. Abstract The aim of the thesis is to investigate the Euler method for GPUs. The goal is to analyze the memory hierarchies, local data-reuse, memory bandwidth of the GPU and the synchronization of the host and the device. Two implemenations are considered closer the diamond tiling and the linear methode, since they are well suited to investigate further improvements. It turns out that the optimizations as local data-reuse and optimum access to the memory bandwidth behave contrary. In the end a combined system (the honeycomb tiling) that combines the advantages of the linear method and the diamond tiling leads to the best results.