Implementierung eines GPU-beschleunigten Kalman-Filters ... · Implementierung eines...

Implementierung einesGPU-beschleunigten

Kalman-Filters mittels OpenCL

Track Fitting für das ATLAS-Experiment am CERN

Masterarbeit

vorgelegt von: Maik Dankel, B.Sc.

Matrikelnummer: 58 28 28

Erstgutachter: Prof. Dr. rer. nat. Nikolaus Wulff

Zweitgutachter: Dr. Sebastian Fleischmann

Datum: 15.10.2013

c© 2013

Implementierung eines GPU-beschleunigten Kalman-Filtersmittels OpenCL

Dieses Werk einschließlich seiner Teile ist urheberrechtlich geschützt. Jede Ver-wertung ausserhalb der engen Grenzen des Urheberrechtgesetzes ist ohne Zustim-mung des Autors unzulässig und strafbar. Das gilt insbesondere für Vervielfältigun-gen, Übersetzungen, Mikroverfilmungen sowie die Einspeicherung und Verarbeitungin elektronischen Systemen.

c© Maik Dankel


Zusammenfassung

Diese Masterarbeit befasst sich mit der Frage, ob die Portierung eines Algo-rithmus zur Spurrekonstruktion auf GPUs einen Performancegewinn gegenüberkonventionellen CPU-Algorithmen erzielen kann. Es werden verschiedene Kon-zepte und der technologische Hintergrund erörtert, der nötig ist, um das Pro-grammiermodell von GPGPU Anwendungen zu verstehen. Anschließend wirdein Kalman-Filter mittels der, von der Khronos Group entwickelten, OpenCLAPI implementiert. Es werden Zeitmessungen vorgenommen, die die Laufzeitdes Algorithmus mit denen einer CPU-Implementierung vergleicht. Die Ergeb-nisse zeigen, dass der Algorithmus, der sich der GPGPU Technologie bedient,lediglich 37% der ursprünglichen Ausführungszeit benötigt.

Abstract

This Master’s thesis deals with the question if porting a track-fitting algo-rithm to GPUs can obtain a better performance than running a conventionalCPU algorithm. First various concepts and the technological background ne-cessary to understand the programming model are discussed in detail. Then aKalman-Filter is implemented using the OpenCL API introduced by KhronosGroup. The execution time of the algorithm is measured and compared withthe execution time of a CPU algorithm. The results show that the GPGPUalgorithm only takes 37% of the original calculation time.

c© Maik Dankel


Danksagung

An dieser Stelle möchte ich mich herzlich bei all denen bedanken, die mich beider Anfertigung dieser Masterarbeit unterstützt haben.

Der größte Dank gilt Herrn Prof. Dr. rer. nat. Nikolaus Wulff, der es mir ermög-licht hat diese Arbeit anzufertigen und der den Kontakt zu der WuppertalerGruppe des ATLAS-Experiments hergestellt hat.

Weiter bin ich Herrn Dr. Sebastian Fleischmann für die durchgehend freundli-che und geduldige Beratung und Betreuung während des Projektes zu großemDank verpflichtet. An dieser Stelle sei auch der restlichen ATLAS-Forschungs-gruppe aus Wuppertal, insbesondere Herrn Prof. Dr. Peter Mättig, für diefreundliche Zusammenarbeit gedankt.

Auch meinem Projektpartner Rene Böing danke ich für die konstruktive Zu-sammenarbeit und auch das Korrekturlesen meiner Arbeit. An dieser Stelle seiauch dem Masterprojekt-Team, bestehend aus Philipp Schoppe und MatthiasTöppe, für ihre Unterstützung gedankt.

Auch meinen weiteren Korrekturlesern bin ich für ihre Hilfe sehr dankbar.Außerdem allen Weiteren, die direkt oder indirekt bei der Erstellung dieserArbeit beteiligt waren.

c© Maik Dankel


Inhaltsverzeichnis

Inhaltsverzeichnis

Abbildungsverzeichnis III

Tabellenverzeichnis IV

Verzeichnis der Listings V

1. Einleitung 11.1. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2. Ziel der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2. GPGPU 42.1. Graphic Processing Units . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1. Architektur . . . . . . . . . . . . . . . . . . . . . . . . . 52.1.2. SMX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.2. OpenCL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2.1. Host-Code . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.2. Device-Code / Kernel . . . . . . . . . . . . . . . . . . . . 182.2.3. Speichermodell . . . . . . . . . . . . . . . . . . . . . . . 212.2.4. Asynchrone Verarbeitung und Synchronisation . . . . . . 23

3. Kalman-Filter 273.1. Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.2. Startwerte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.3. Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4. Implementierung 334.1. Data-IO und Datenstrukturen . . . . . . . . . . . . . . . . . . . 334.2. Host-Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.3. Device-Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.3.1. Kalman-Filter . . . . . . . . . . . . . . . . . . . . . . . . 404.3.2. Smoothing . . . . . . . . . . . . . . . . . . . . . . . . . . 53

c© Maik Dankel I


Inhaltsverzeichnis4.4. Optimierungen . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.4.1. Track-Parallelisierung . . . . . . . . . . . . . . . . . . . . 554.4.2. Asynchrone Verarbeitung . . . . . . . . . . . . . . . . . . 574.4.3. Host-Parallelisierung . . . . . . . . . . . . . . . . . . . . 604.4.4. Numerische Stabilität und Symmetrie . . . . . . . . . . . 604.4.5. Matrixinvertierung . . . . . . . . . . . . . . . . . . . . . 61

5. Performance-Vergleich 625.1. Optimierungsstufen . . . . . . . . . . . . . . . . . . . . . . . . . 625.2. OpenCL / CUDA / CPU . . . . . . . . . . . . . . . . . . . . . . 63

6. Fazit 65

7. Ausblick 66

Literaturverzeichnis 67

Eidesstattliche Erklärung 69

A. Anhang i

c© Maik Dankel II


Abkürzungsverzeichnis

Abbildungsverzeichnis

1.1. Schematische Darstellung des ATLAS Detektors. . . . . . . . . . 2

2.1. Komplettes Blockdiagramm des GK110. . . . . . . . . . . . . . 52.2. Blockdiagramm einer SMX-Einheit eines GK110-Chips. . . . . . 72.3. Arbeitsweise des Warp Schedulers . . . . . . . . . . . . . . . . . 82.4. Aufbau einer OpenCL Plattform . . . . . . . . . . . . . . . . . . 92.5. Veranschaulichung von Global- und Local-ID . . . . . . . . . . . 192.6. Veranschaulichung von Work-Groups, Work-Items und Dimen-

sionen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212.7. OpenCL Speichermodell . . . . . . . . . . . . . . . . . . . . . . 222.8. Synchronisierung von Work-Items . . . . . . . . . . . . . . . . . 24

3.1. Der rekursive Ablauf des Kalman-Filters. Der time update Schrittbildet den aktuellen Zustand auf einen geschätzten zukünfti-gen Wert ab. Der measurement update Schritt korrigiert dengeschätzten Wert mittels einer aktuellen Messung. . . . . . . . . 28

3.2. Schematische Darstellung von Teilchenspur und verfälschter Mes-sung innerhalb eines Detektors . . . . . . . . . . . . . . . . . . . 29

3.3. Schematische Darstellung der durch den Kalman-Filter korri-gierten Spur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.1. Arbeitsweise der Implementation . . . . . . . . . . . . . . . . . 364.2. Speicheraufteilung . . . . . . . . . . . . . . . . . . . . . . . . . 374.3. Track Parallelisierung . . . . . . . . . . . . . . . . . . . . . . . . 554.4. Speicheraufteilung . . . . . . . . . . . . . . . . . . . . . . . . . 56

c© Maik Dankel III


Tabellenverzeichnis

Tabellenverzeichnis

4.1. Hilfsfunktionen Matrix-/Vektor-Multiplikation . . . . . . . . . . 44

5.1. Spezifikationen des Testsystems . . . . . . . . . . . . . . . . . . 625.2. Eigenschaften des verwendeten Testdatensatzes . . . . . . . . . 625.3. Performancevergleich Optimierungsstufen . . . . . . . . . . . . . 635.4. Ausführungszeiten der Implementierungen . . . . . . . . . . . . 645.5. Performancevergleich der Implementierungen mit Benchmarks . 645.6. Performancevergleich der Implementierungen mit Benchmark 500 64

c© Maik Dankel IV


Verzeichnis der Listings

Verzeichnis der Listings

2.1. oclGetPlatformID() Signatur . . . . . . . . . . . . . . . . . . . . 102.2. clGetDeviceIDs() Signatur . . . . . . . . . . . . . . . . . . . . . 102.3. clCreateContext() Signatur . . . . . . . . . . . . . . . . . . . . . 112.4. clCreateCommandQueue() Signatur . . . . . . . . . . . . . . . . 122.5. clCreateBuffer() Signatur . . . . . . . . . . . . . . . . . . . . . . 132.6. clEnqueueWriteBuffer() Signatur . . . . . . . . . . . . . . . . . 142.7. clCreateProgramWithSource() Signatur . . . . . . . . . . . . . . 152.8. clBuildProgram() Signatur . . . . . . . . . . . . . . . . . . . . . 152.9. clCreateKernel() Signatur . . . . . . . . . . . . . . . . . . . . . 162.10. clSetKernelArg() Signatur . . . . . . . . . . . . . . . . . . . . . 162.11. clEnqueueNDRangeKernel() Signatur . . . . . . . . . . . . . . . 162.12. clEnqueueReadBuffer() Signatur . . . . . . . . . . . . . . . . . . 172.13. OpenCL clean-up Methoden . . . . . . . . . . . . . . . . . . . . 182.14. OpenCL Synchronisationspunkt barrier() . . . . . . . . . . . . . 242.15. clWaitForEvents() Signatur . . . . . . . . . . . . . . . . . . . . 252.16. clFinish() Signatur . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.1. Definition eines Events . . . . . . . . . . . . . . . . . . . . . . . 344.2. Definition eines Tracks . . . . . . . . . . . . . . . . . . . . . . . 344.3. Definition eines Tracks . . . . . . . . . . . . . . . . . . . . . . . 344.4. Pseudocode Filtering-Prozess . . . . . . . . . . . . . . . . . . . 374.5. Überprüfung der Dimensionalität . . . . . . . . . . . . . . . . . 384.6. Invertierung der Kovarianzmatrizen . . . . . . . . . . . . . . . . 394.7. Pseudo-Code Filtering-Kernel . . . . . . . . . . . . . . . . . . . 404.8. Funktionskopf des doFiltering Kernels . . . . . . . . . . . . . . . 414.9. Start des Vorwärts-Fits . . . . . . . . . . . . . . . . . . . . . . . 424.10. Time-Update Hilfsfunktion . . . . . . . . . . . . . . . . . . . . . 434.11. Hilfsfunktion zur Matrix-Vektor-Multiplikation . . . . . . . . . . 444.12. Erster Teil der Berechnung von Ck|k−1 . . . . . . . . . . . . . . 454.13. Zweiter Teil der Berechnung von Ck|k−1 . . . . . . . . . . . . . . 45

c© Maik Dankel V


Verzeichnis der Listings4.14. Berechnung von Kk im eindimensionalen Fall . . . . . . . . . . . 474.15. Kalman-Gain: Invertierung 2D . . . . . . . . . . . . . . . . . . . 484.16. Kalman-Gain: Matrixmultiplikation 2D . . . . . . . . . . . . . . 494.17. Berechnung pk|k 1D . . . . . . . . . . . . . . . . . . . . . . . . . 504.18. Berechnung pk|k 2D . . . . . . . . . . . . . . . . . . . . . . . . . 504.19. Berechnung Ck|k 1D . . . . . . . . . . . . . . . . . . . . . . . . . 514.20. Berechnung Ck|k 1D . . . . . . . . . . . . . . . . . . . . . . . . . 524.21. Smoothing Kernel . . . . . . . . . . . . . . . . . . . . . . . . . . 544.22. Verwendung von Work-Groups . . . . . . . . . . . . . . . . . . . 574.23. Datenstruktur für die Nutzung des OpenCL Eventmodells . . . 584.24. Pseudocode asynchrone Verarbeitung . . . . . . . . . . . . . . . 59

A.1. Minimalbeispiel einer OpenCL Implementierung . . . . . . . . . i

c© Maik Dankel VI


1. Einleitung

1. Einleitung

Das europäische Kernforschungszentrum CERN (Conseil Européen pour la Re-cherche Nucléaire) mit Sitz in Meyrin (Kanton Genf) ist eine Großforschungs-einrichtung, die sich mit physikalischer Grundlagenforschung beschäftigt. Ins-besondere werden mit Hilfe des "Großen Hadronen-Speicherrings" (engl.: Lar-ge Hadron Collider, kurz: LHC), dem weltgrößten Teilchenbeschleuniger, diegrundlegenden Bestandteile der Materie erforscht. Dort werden Teilchen mitannähernd Lichtgeschwindigkeit zur Kollision gebracht, was den PhysikernAufschluss darüber gibt, wie diese Teilchen miteinander wechselwirken. Au-ßerdem erhofft man sich neue Erkenntnisse über die fundamentalen Naturge-setze[CERN].Dem LHC angeschlossen sind sieben Experimente, die Teilchendetektoren be-treiben, um die unzähligen, bei der Kollision erzeugten Teilchen, zu messenund zu analysieren. Diese werden von Kollaborationen von Wissenschaftlernweltweit betrieben. Die beiden größten Experimente sind ATLAS und CMS,die unabhängig voneinander konzipiert wurden, was unverzichtbar ist, um neueEntdeckungen validieren zu können.

1.1. Motivation

ATLAS ist eines der großen Experimente am LHC des CERN. Der Detektorist mit einer Länge von etwa 45 Metern und einem Durchmesser von 25 Meternder größte Detektor am LHC (siehe Abbildung 1.1). Ziel des Experiments istdie Suche nach unentdeckten Teilchen, wie zum Beispiel dem Higgs-Boson, so-wie die Untersuchung einer womöglich vorhandenen Substruktur von Leptonenund Quarks, beides Klassen von Elementarteilchen aus dem Standardmodellder Teilchenphysik.Um diese Teilchen zu finden und zu analysieren werden Protonen oder Blei-Ionen zu Strahlen gebündelt und bei nahezu Lichtgeschwindigkeit zur Kollisi-on gebracht, um die Wechselwirkung dieser Teilchen bei sehr hohen Energien

c© Maik Dankel 1


1. Einleitung

Abbildung 1.1.: Schematische Darstellung des ATLAS Detektors.Quelle: [ATLAS]

(nach dem Anfang 2013 eingeleiteten, etwa 2-jährigen technischen Stopp (LS1),bis zu 14TeV Schwerpunktenergie) zu untersuchen.Bei der Kollision entsteht eine Wechselwirkung zwischen den kollidierten Teil-chen und diese zerfallen in weitere Subteilchen. Die einzelnen Sekundärteilchenzerfallen ihrerseits wieder in weitere Teilchen. Die so entstehenden Teilchen-schauer erzeugen Myriaden von Messpunkten, die es auszuwerten gilt.Das ATLAS-Experiment erzeugt mit den Messungen dieser Teilchenschauergewaltige Datenmengen von über 300 MB pro Sekunde [ATLAS], die auszu-werten aktuell sehr viel Zeit in Anspruch nimmt. Aus diesem Grund wirdmomentan viel in die Forschung und Entwicklung neuer, effizienterer Wegeinvestiert diese Daten zu analysieren. Ein Ansatz ist es, die Berechnungen (zu-mindest teilweise) zu parallelisieren und, statt wie bisher auf CPUs, auf GPUsauszuführen. Diese sind in der Lage mehrere tausend Berechnungen zeitgleich(echt parallel) durchzuführen, was bei bestimmten Aufgaben einen erheblichenGeschwindigkeitsgewinn bewirken kann.

1.2. Ziel der Arbeit

Ziel dieser Masterarbeit ist es, einen Algorithmus zur Spurrekonstruktion (indiesem Fall ein Kalman-Filter), wie er in Teilen der ATLAS Software eingesetzt

c© Maik Dankel 2


1. Einleitungwird, so zu implementieren, dass er auf GPUs ausgeführt werden kann. Es giltherauszufinden, ob dieses spezielle Problem effizient auf GPUs berechnet wer-den kann. Die Software wird unter Benutzung der freien OpenCL Program-mierschnittstelle implementiert und kontinuierlich hinsichtlich Geschwindig-keit optimiert. Um den Geschwindigkeitszuwachs feststellen zu können, werdenwiederholt Zeitmessungen der Berechnungsdauer durchgeführt. Parallel dazuwird eine CPU-Anwendung entwickelt, die die selben Daten auswertet1 umeinen Referenzwert für die Berechnungsdauer zu erhalten und die Ergebnisseevaluieren zu können.

1Die Implementierung des CPU-Codes war Teil eines Masterprojektes und nicht Bestandteildieser Masterarbeit

c© Maik Dankel 3


2. GPGPU

2. GPGPU

Im Oktober 2010 fand im High-Performance Computing eine Revolution statt.China stellte seinen neuen Supercomputer Tianhe-1A vor, der als erster Com-puter seiner Klasse GPUs als Koprozessoren verwendete. Mit einer Rechen-leistung von 2,566 PFLOPS (1015 floating-point operations per second) ist erfast 50% schneller als der bis dahin schnellste Supercomputer Jaguar (1,759PFLOPS). Eine weitere Besonderheit ist, dass er zwar fast doppelt so rechen-stark ist, jedoch nur knapp 58% der elektrischen Leistung von Jaguar benötigtund damit bedeutend energieeffizienter arbeitet (Vgl. [Scarpino 2012], S. 3f).Der Einsatz von GPUs als Koprozessoren für nicht-graphische Berechnungenwird als general purpose GPU computing oder GPGPU computing bezeichnetund ist im Bereich High-Performance Computing momentan immer mehr imVormarsch.

2.1. Graphic Processing Units

Grafikprozessoren können in speziellen Bereichen erheblich mehr Leistung brin-gen als herkömmliche CPUs. Dies gilt jedoch nur unter bestimmten Vorausset-zungen. Um nachzuvollziehen unter welchen Umständen sie leistungsfähigersein können, muss zunächst die Architektur einer GPU genauer betrachtetwerden. Eine der aktuell am weitesten fortgeschrittenen GPUs wurde 2012von NVIDIA unter dem Namen GK110 vorgestellt. Dieser Chip wurde haupt-sächlich für die NVIDIA Tesla Koprozessoren entwickelt und basiert auf derneuen Kepler Architektur. Ziel der Entwicklung dieses Chips ist es die höchsteParallel-Computing Performance der Welt zu erreichen [NVIDIA]. Eine kom-plette Kepler GK110 Implementierung besteht aus 15 Streaming Multiprozes-soren (SMX) und sechs 64-Bit Memory Controllern. Inzwischen wird dieserChip teilweise in NVIDIA’s Consumer-Modellen, der GeForce Reihe, verbaut.Sowohl bei den GeForce Modellen, als auch den anderen Baureihen, sind je-doch oftmals einige der Komponenten deaktiviert (es steht beispielsweise eine

c© Maik Dankel 4


2. GPGPU

Abbildung 2.1.: Komplettes Blockdiagramm des GK110.Quelle: [NVIDIA] S.6

geringere Anzahl an SMX zur Verfügung).

2.1.1. Architektur

Das in Abbildung 2.1 abgebildete Blockdiagramm zeigt den schematischenAufbau eines GK110 Chips. Nachfolgend werden die einzelnen Elemente zurBegriffsklärung und zum besseren Verständnis kurz erläutert.

• Gigathread EngineDie Gigathread Engine ist dafür zuständig, die zu berechnenden Datenin Threadblöcke aufzuteilen und an die SMX Einheiten weiterzuleiten.Des Weiteren ist sie verantwortlich für Kontextwechsel zwischen den ein-zelnen Blöcken.

• Streaming Multiprozessor (SMX)Die Next Generation Streaming Multiprocessors sind das Herzstück derKepler GPUs. Ein SMX ist vom Prinzip her vergleichbar mit einer CPU

c© Maik Dankel 5


2. GPGPUund ist für die Aufteilung der Blöcke auf mehrere Threads verantwortlich,die auf die einzelnen Kerne verteilt werden. Die einzelnen Komponenteneines solchen Multiprozessesors werden in Unterabschnitt 2.1.2 genauerbehandelt.

• PCI Express 3.0 Host InterfaceÜber die PCIe Schnittstelle ist die Grafikkarte mit dem Hostsystem ver-bunden. PCIe 3.0 x16 erlaubt Datenübertragungsraten von bis zu 16GBpro Sekunde.

• L2 CacheDer Level 2 Cache dient dazu, Daten zwischen den einzelnen Prozessoren(hier den SMX) auszutauschen. Er dient demnach der Verwaltung desglobalen Speichers, der unter den SMX geteilt wird.

• Memory ControllerDer GK110 Chip verfügt über insgesamt sechs 64-Bit Memory Controller,die die Speicherverwaltung übernehmen. Die größten Karten, die denGK110 Chip einsetzen, verfügen über 6GB GDDR5 Speicher.

2.1.2. SMX

Abbildung 2.2 zeigt den schematischen Aufbau eines SMX des GK110-Grafik-chips. Auch hier werden die einzelnen Komponenten nachfolgend erläutert undbeschrieben.

• KerneDer große mittlere Teil in Abbildung 2.2 stellt die Rechenkerne des Multi-prozessors dar. Jede SMX-Einheit besteht aus 192 Single-Precision- und64 Double-Precision Kernen, die für mathematische Berechnungen ge-nutzt werden können. Außerdem existieren je 32 Load/Store Einheitenund sogenannte Special Function Units, die der Berechnung speziellermathematischer Funktionen wie Sinus, Cosinus, Wurzeln, etc. dienen.

• Warp SchedulerEin SMX teilt Threads in Gruppen von 32 Threads ein, die Warps ge-nannt werden. Jede der SMX-Einheiten besitzt vier Warp Scheduler, diees erlauben vier Warps gleichzeitig zu starten und auszuführen. Außer-dem sind pro Warp Scheduler zwei Dispatch Einheiten vorhanden, die

c© Maik Dankel 6


2. GPGPU

Abbildung 2.2.: Blockdiagramm einer SMX-Einheit eines GK110-Chips.Quelle: [NVIDIA] S.8

es ermöglichen zwei unabhängige Befehle pro Takt an einen Warp zuschicken. Abbildung 2.3 veranschaulicht dieses Verhalten.

• Instruction CacheDer Instruction Cache ist ein spezieller Speicher zum Zwischenspeichernvon Befehlen, bevor diese an die einzelnen Warps verteilt werden.

• Shared Memory / L1-CacheJede SMX-Einheit verfügt über 64 KB lokalen Speicher, der sowohl alsshared Memory innerhalb des SMX dienen kann, als auch als Level 1-Cache. Dieser Speicher kann entweder als 16 KB L1-Cache und 48 KB

c© Maik Dankel 7


2. GPGPU

Abbildung 2.3.: Arbeitsweise des Warp SchedulersQuelle: [NVIDIA] S.10

Shared Memory, 48 KB L1-Cache und 16 KB Shared Memory, odergleichmäßig auf je 32 KB aufgeteilt werden.

• Read-Only Data CacheDer Read-Only Data Cache ist ein spezieller Speicher für Daten, auf dienur lesend zugegriffen werden darf. Die Speichergröße beträgt 48 KB.

2.2. OpenCL

OpenCL ist eine Programmierschnittstelle (API) zum Schreiben von Softwa-re, die plattformunabhängig ausgeführt werden kann. Der OpenCL-Code kannlaut Spezifikation auf CPUs, GPUs, DSPs und anderen Prozessoren ausgeführtwerden. OpenCL wurde zunächst von Apple in Kooperation mit führenden Fir-men aus der Industrie2 entwickelt und Mitte 2008 zur Standardisierung bei derKhronos Group eingereicht. Ende 2008 wurde dann die OpenCL 1.0 Spezifika-tion veröffentlicht [KHRONOS].

OpenCL stellt eine Top-Level-Abstraktion für Low-Level Hardware Routinendar, um hochparallele Anwendungen zu implementieren. Dabei kann sowohl

2 z.B. 3DLABS, Activision Blizzard, AMD, ARM, Electronic Arts, Ericsson, IBM, IntelCorporation, Motorola, Nokia, NVIDIA, Samsung, Texas Instruments

c© Maik Dankel 8


2. GPGPU

Abbildung 2.4.: Aufbau einer OpenCL PlattformQuelle: [KHRONOS 2008] S.13

Datenparallelität, als auch Task-Parallelität erreicht werden. Der eigentlicheOpenCL-Code wird in OpenCL C implementiert, eine auf ISO C99 basierendeSprache, die um einige Datentypen und Funktionen zur parallelen Verarbeitungerweitert wurde. Allerdings ist OpenCL C an einigen Stellen eingeschränkt ge-genüber dem C99-Standard. So existieren beispielsweise keine Funktionszeiger,Arrays müssen eine feste Größe besitzen und Rekursion ist nicht möglich.

OpenCL Software ist in zwei Komponenten aufgeteilt, zum einen der Host-Code und auf der anderen Seite der Device-Code. Als Host wird im Folgendenein System bezeichnet, das ein oder mehrere OpenCL-fähige Geräte (Devi-ces) enthält und von dem/denen aus die Software ausgeführt wird. Ein Devicewiederum ist eine Grafikkarte oder ein Koprozessor, auf dem der eigentlicheOpenCL-Code ausgeführt wird.Im grundsätzlichen Aufbau einer OpenCL Plattform verfügt ein Host über einoder mehrere Devices und steuert diese an (siehe Abbildung 2.4). Jedes Deviceverfügt über mehrere Recheneinheiten (siehe Unterabschnitt 2.1.1, SMX), wel-che wiederum über mehrere Rechenkerne verfügen (siehe Unterabschnitt 2.1.2,Kerne).

In den folgenden Abschnitten wird erläutert, welche Schritte nötig sind umlauffähigen OpenCL-Code zu erzeugen und einige Besonderheiten dabei werdenhervorgehoben.

c© Maik Dankel 9


2. GPGPU2.2.1. Host-Code

Hostseitig werden sämtliche, für die OpenCL-Laufzeitumgebung wichtigen,Komponenten implementiert und angesteuert. Der Hostcode als solches kann inregulärem C oder C++ geschrieben werden, allerdings gibt es einige OpenCL-spezifische Datenstrukturen und Funktionsaufrufe, auf die nachfolgend genauereingegangen wird.

Erstellen der grundlegenden OpenCL-Umgebung

Dazu zählt als Erstes die sogenannte Plattform, welche durch eine cl_platform_id repräsentiert wird. Eine Plattform ist ein Host mit 1-n Devices (siehe Ab-bildung 2.4). Mit folgendem Aufruf wird ein cl_platform Objekt erzeugt:

1 // Returns an error code2 cl_int oclGetPlatformID (cl_platform_id ∗platform) // Pointer to the platform object

Listing 2.1: oclGetPlatformID() Signatur

Es wird ein Zeiger auf eine cl_platform_id übergeben, welcher durch den Auf-ruf gesetzt wird und damit die Referenz zur Plattform herstellt. Als Rückgabe-wert wird ein OpenCL-weit eindeutiger numerischer Wert zurückgegeben, dereinen Fehlercode darstellt.

Sobald Zugriff auf die Plattform besteht, kann auf alle Devices zugegriffenwerden, die mit dieser Plattform assoziiert sind. Dazu muss die Referenz zu ei-nem oder mehreren Devices hergestellt werden, was über den folgenden Aufrufgeschieht:

1 cl_int clGetDeviceIDs (2 cl_platform_id platform,3 cl_device_type device_type, // bitfield for choosing the device type4 cl_uint num_entries, // The number of cl_device entries that can be added to

\emph{devices}5 cl_device_id ∗devices, // Pointer to the device object6 cl_uint ∗num_devices) // number of devices matching the device_type

Listing 2.2: clGetDeviceIDs() Signatur

c© Maik Dankel 10


2. GPGPUHier wird festgelegt wie viele und was für Devices für die angegebene Plattformangesprochen werden. Mögliche Typen für device_type sind:

• CL_DEVICE_TYPE_GPU

• CL_DEVICE_TYPE_CPU (Host-Prozessor, Single- oder Multicore)

• CL_DEVICE_TYPE_ACCELERATOR (dedizierter OpenCL Kopro-zessor)

• CL_DEVICE_TYPE_DEFAULT

• CL_DEVICE_TYPE_ALL

Nach einem erfolgreichen Aufruf sind die Devices eindeutig über deren cl_device_id definiert und können darüber angesprochen werden. Ebenso ist dann dieAnzahl der Komponenten, die dem angegebenen Typ entsprechen, in num _de-vices gespeichert.

Nach Definition einer Plattform und ein oder mehrerer Devices, wird ein Con-text erzeugt. Dieser definiert die gesamte OpenCL Umgebung inklusive derOpenCL-Devices, Kernel, Command-Queues, etc. und dient dazu OpenCL-Code auf den jeweiligen Devices auszuführen. Ein Context wird über den nach-folgenden Funktionsaufruf erzeugt:

1 cl_context clCreateContext (2 const cl_context_properties ∗properties, // properties bitfield (see

specification )3 cl_uint num_devices, // Number of devices4 const cl_device_id ∗devices, // Pointer to the devices object5 void (∗pfn_notify)(const char ∗errinfo , const void ∗private_info, size_t cb,

void ∗user_data), // callback function to retrieve errors6 void ∗user_data, // Passed as user_data argument when pfn_notify is called7 cl_int ∗errcode_ret) // error code result

Listing 2.3: clCreateContext() Signatur

Die wichtigen Parameter dieses Funktionsaufrufs sind num_devices (wie vieleDevices werden mit diesem Context assoziiert) und *devices, ein Zeiger, derdie genauen Devices angibt und durch den vorherigen Aufruf von clDetDevi-ceIDs() gesetzt wird. Zusätzlich kann eine callback-Funktion angegeben wer-den, die von der OpenCL-Implementierung genutzt wird, um Informationen

c© Maik Dankel 11


2. GPGPUüber Fehler zu liefern, die in diesem Context auftreten. Soll diese Funktio-nalität nicht genutzt werden, können die Argumente (*pfn_notify)(...) und*user_data NULL gesetzt werden3.

Abschließend wird eine Command-Queue erzeugt. In dieser Warteschlan-ge werden Anweisungen, die auf dem jeweiligen Device ausgeführt werdensollen, abgelegt und von diesem letztendlich abgearbeitet. Erzeugt wird eineCommand-Queue mit dem Aufruf von clCreateCommandQueue():

1 cl_command_queue clCreateCommandQueue (2 cl_context context,3 cl_device_id device,4 cl_command_queue_properties properties, // bitfiled for properties5 cl_int ∗errcode_ret) // error code result

Listing 2.4: clCreateCommandQueue() Signatur

Zwischen Devices und Command-Queues besteht eine 1:n Beziehung. EinerWarteschlange ist über das Argument device immer genau ein OpenCL-Devicezugeordnet, jedoch kann ein Device mehrere Queues haben, was unter Um-ständen einen höheren Grad an Parallelisierung ermöglichen kann (Vgl. dazuUnterabschnitt 2.2.4).Die Anweisungen, die in der Queue abgelegt werden, können entweder ”in-order” oder ”out-of-order” ausgeführt werden. Dies wird über das propertiesArgument festgelegt. Standardmäßig werden die Befehle in der Reihenfolge ab-gearbeitet, in der sie dort abgelegt werden (first in, first out). Soll dieses Verhal-ten geändert werden, kann für properties CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE gesetzt werden. Ist das in-order Verhalten ge-wünscht, kann hierfür 0 eingesetzt werden.

Speicher allozieren

Mit den vorherigen Aufrufen ist der Host grundlegend konfiguriert. Es könnenDevices angesprochen und Anweisungen an die Command-Queue(s) weiterge-geben werden. Als Nächstes wird Speicherplatz auf dem Device zur Verfügunggestellt. Das erfolgt mit folgendem Funktionsaufruf:

3selbiges gilt auch für *properties

c© Maik Dankel 12


2. GPGPU

1 cl_mem clCreateBuffer (2 cl_context context, // The context where the memory will be allocated3 cl_mem_flags flags,4 size_t size , // The size in bytes5 void ∗host_ptr,6 cl_int ∗errcode_ret)

Listing 2.5: clCreateBuffer() Signatur

Dieser Aufruf erzeugt ein cl_mem Objekt, das einen Speicherbereich repräsen-tiert. Dieser kann verschiedene Eigenschaften haben, je nachdem welche Wertefür das Argument flags gesetzt werden. Mögliche Flags sind

• CL_MEM_READ_WRITE

• CL_MEM_WRITE_ONLY

• CL_MEM_READ_ONLY

• CL_MEM_COPY_HOST_PTR

• CL_MEM_ALLOC_HOST_PTR

• CL_MEM_USE_HOST_PTR

Mit den ersten drei Flags können Zugriffsrechte für den Speicher angegebenwerden. Nur lesender, nur schreibender Zugriff, oder beides. Wird keines ge-setzt, so kann standardmäßig sowohl lesend, als auch schreibend zugegriffenwerden. Die letzten drei Flags geben an, ob und wie Daten (beim Erstellen desSpeichers) zwischen Host und Device übertragen werden. Den einfachsten Fallstellt CL_MEM_COPY_ HOST_PTR dar, da hier Speicher auf dem Devi-ce alloziert und gleichzeitig mit den Daten, die durch host_pointer adressiertsind, initialisiert wird. Mit CL_MEM_USE_HOST_PTR wird Speicher ineinem Bereich alloziert, der sowohl vom Host, als auch vom Device erreichbarist. Hierbei ist Vorsicht geboten, da der Speicher hostseitig manipuliert wird,sobald der Kernel gestartet wird. Das Flag CL_MEM_ ALLOC_HOST_PTRbewirkt, dass hostseitig sogenannter Pinned-Memory erzeugt wird. Dies ist einSpeicherbereich, der nicht auf Sekundärspeicher (z.B. die Festplatte) ausgela-gert wird. Das kann erheblich höhere Datentransferraten ermöglichen. Wirdkeines der letzten drei Flags angegeben, so wird lediglich Speicher auf demDevice reserviert. Jeweils eines von den ersten und letzten drei Flags kannbitweise verodert werden.

c© Maik Dankel 13


2. GPGPUSofern nicht das Flag CL_MEM_USE_HOST_PTR (bei dem Hostspeichergenutzt wird), bzw. CL_MEM_COPY_HOST_PTR (wobei implizit Datenauf das Device kopiert werden) gesetzt wird, muss, sofern Daten zum Devicetransferiert werden sollen, ein expliziter Kopiervorgang angestoßen werden.Das geschieht mit clEnqueueWriteBuffer():

1 cl_int clEnqueueWriteBuffer (2 cl_command_queue command_queue,3 cl_mem buffer, // to which buffer4 cl_bool blocking_write, // determines if this call is blocking or not5 size_t offset , // offset from the beginning6 size_t cb, // size to be written (in bytes)7 const void ∗ptr, // pointer to host memory8 cl_uint num_events_in_wait_list,9 const cl_event ∗event_wait_list,

10 cl_event ∗event)

Listing 2.6: clEnqueueWriteBuffer() Signatur

Hierbei wird angegeben welcher Speicherbereich (bestimmt durch *ptr und cb)an welche Stelle auf dem Device (angegeben durch das vorher erzeugte cl_memObjekt, hier buffer) kopiert wird. Wie der Name der Funktion bereits impli-ziert, wird das Kommando in die Command-Queue gehängt und kann über dasboolsche Argument blocking_write blockierend oder nicht blockierend ausge-führt werden. Sofern blocking_write auf FALSE gesetzt ist, muss mittels desletzten Arguments event sichergestellt werden, dass der Aufruf beendet ist,bevor auf diesen Daten gearbeitet werden kann. Sowohl clEnqueueWriteBuf-fer, als auch alle noch folgenden clEnqueue* Befehle arbeiten nach diesemPrinzip.4

Programm und Kernel

Sowohl Programm und Kernel enthalten ausführbaren Code. Ein Kernel (siehedazu Unterabschnitt 2.2.2) ist allerdings eine einzelne Funktion, die auf einemDevice ausgeführt wird. Im Gegensatz dazu ist ein Programm eine Sammlungvon mehreren Kerneln.

4mehr zum Thema Events in Unterabschnitt 2.2.4

c© Maik Dankel 14


2. GPGPUEin Programm wird in OpenCL durch eine cl_program Datenstruktur reprä-sentiert, die mit folgendem Funktionsaufruf erzeugt wird:

1 cl_program clCreateProgramWithSource (2 cl_context context,3 cl_uint count, // number of files4 const char ∗∗strings , // array of strings , each one is a file5 const size_t ∗lengths, // array specifying the file lengths6 cl_int ∗errcode_ret) // error code to be returned

Listing 2.7: clCreateProgramWithSource() Signatur

Dieser Aufruf erzeugt ein Programm aus einer oder mehreren Dateien (angege-ben durch das Argument count) und erwartet ein Array von Strings, das denInhalt dieser Dateien enthält. Das heißt, eine Datei muss zunächst ausgelesenund zwischengespeichert werden. Diese Datei enthält den eigentlichen Kernel-Code.Ist dies erfolgt, wird das Programm kompiliert. Das erfolgt mit dem Aufrufvon clBuildProgram(). Die Signatur sieht folgendermaßen aus:

1 cl_int clBuildProgram (2 cl_program program,3 cl_uint num_devices,4 const cl_device_id ∗device_list,5 const char ∗options, // Compiler options, see the specifications for more

details6 void (∗pfn_notify)(cl_program, void ∗user_data),7 void ∗user_data)

Listing 2.8: clBuildProgram() Signatur

Diese Funktion benötigt sowohl das zuvor erzeugte cl_program, als auch dieAnzahl der Devices und die Referenz zu diesen. Mit dem Argument optionskönnen Compileroptionen angegeben werden. Ähnlich wie beim Aufruf vonclCreateContext() existiert auch hier eine Callback-Funktion. Sofern diese an-gegeben wird, wird sie gerufen, sobald der Build-Vorgang abgeschlossen ist.Außerdem blockiert der Aufruf von clBuildProgram() in diesem Fall nicht undes wird mit der Ausführung des Codes direkt fortgefahren. Ist dies nicht ge-wünscht, kann für die beiden letzten Argumente NULL eingesetzt werden.

c© Maik Dankel 15


2. GPGPUAbschließend wird mittels clCreateKernel() ein Kernel erzeugt, der später alsEinstiegspunkt zu dem auf der GPU ausführbaren Code dient.

1 cl_kernel clCreateKernel (2 cl_program program, // The program where the kernel is3 const char ∗kernel_name, // The name of the kernel4 cl_int ∗errcode_ret)

Listing 2.9: clCreateKernel() Signatur

Dieser Aufruf benötigt das vorher erzeugte und kompilierte Programm sowieden Namen des Kernels, der geladen werden soll. Hiermit ist nicht der Da-teiname, sondern der tatsächliche Name der Kernelfunktion gemeint, die imCode aufgerufen wird.

Da ein Kernel gewissermaßen eine Funktion darstellt, benötigt dieser üblicher-weise Argumente, die ihm übergeben werden. Dazu ist ein weiterer Aufrufnötig, der wie folgt aussieht:

1 cl_int clSetKernelArg (2 cl_kernel kernel , // a valid kernel3 cl_uint arg_index, // index of the argument, goes from 0 to n−14 size_t arg_size, // size of the argument value5 const void ∗arg_value) // pointer to data that should be used as the argument

Listing 2.10: clSetKernelArg() Signatur

Dieser muss für n Argumente n-mal aufgerufen werden. Es wird der Kernelangegeben, für den diese Argumente bestimmt sind, sowie ein Index, der dasjeweilige Argument angibt. Bei beispielsweise n=3 Argumenten muss dieseFunktion drei mal aufgerufen werden, mit arg_index ∈ {0, 1, 2}. Die Argumen-te arg_size und arg_value geben die Größe und das zu übergebene Argumentan. Üblicherweise sind das ein cl_mem Objekt, das zuvor mit clCreateBuffer()erstellt wurde und dessen Größe (sizeof(cl_mem)).

Um die Ausführung des Kernels zu starten ist der nachfolgende Aufruf vonclEnqueueNDRangeKernel() nötig.

1

2 cl_int clEnqueueNDRangeKernel (3 cl_command_queue command_queue, // A valid command−queue4 cl_kernel kernel , // A valid kernel

c© Maik Dankel 16


2. GPGPU5 cl_uint work_dim, // Choose if using 1D, 2D or 3D Work−Items and Work−

Groups6 const size_t ∗global_work_offset, // Mmust be NULL, not supported by

OpenCL yet7 const size_t ∗global_work_size, // The total number of Work−Items (must have

work_dim dimensions)8 const size_t ∗local_work_size, // The number of Work−Items per Work−Group

(must have work_dim dimensions)9 cl_uint num_events_in_wait_list,

10 const cl_event ∗event_wait_list,11 cl_event ∗event)

Listing 2.11: clEnqueueNDRangeKernel() Signatur

Es wird eine Command-Queue sowie der auszuführende Kernel übergeben.Außerdem muss die Dimension angegeben werden, in der die Work-Items ar-rangiert sind (zur Erläuterung siehe Unterabschnitt 2.2.2). global_work_sizegibt die Gesamtanzahl der Work-Items an, local_work_size gibt die Anzahlan Work-Items pro Work-Group an und muss ein Teiler von global_work_sizesein. Wird für local_work_size NULL eingesetzt berechnet die OpenCL Im-plementierung selbsständig eine geeignete Aufteilung in Work-Groups. NachAufruf von clEnqueueNDRangeKernel() wird der Befehl zur Ausführung desKernels in die Command-Queue eingehängt und anschließend gestartet.

Speicher auslesen

Nach Ausführen des Kernels werden die Daten vom Devicespeicher ausgelesen.Dazu existiert eine weitere Funktion mit Namen clEnqueueReadBuffer():

1 cl_int clEnqueueReadBuffer (2 cl_command_queue command_queue,3 cl_mem buffer, // from which buffer4 cl_bool blocking_read, // determines if this call is blocking or not5 size_t offset , // offset from the beginning6 size_t cb, // size to be read (in bytes)7 void ∗ptr, // pointer to the host memory8 cl_uint num_events_in_wait_list,9 const cl_event ∗event_wait_list,

10 cl_event ∗event)

Listing 2.12: clEnqueueReadBuffer() Signatur

c© Maik Dankel 17


2. GPGPUDieser Aufruf erfolgt analog zu clEnqueueWriteBuffer(). In diesem Fall gibtbuffer an aus welchem Speicher gelesen wird und *ptr ist eine hostseitige Spei-cheradresse, an der die Daten abgelegt werden.

Clean up

Abschließend werden jegliche zuvor erzeugten Objekte freigegeben. Dies mussexplizit geschehen, da die OpenCL Implementierung sonst dauerhaft Ressour-cen belegt.

1 cl_int clReleaseMemObject ( cl_mem memobj)2 cl_int clReleaseProgram (cl_program program)3 cl_int clReleaseKernel (cl_kernel kernel)4 cl_int clReleaseCommandQueue(cl_command_queue command_queue)5 cl_int clReleaseContext (cl_context context)

Listing 2.13: OpenCL clean-up Methoden

Damit sind alle wichtigen Aufrufe erläutert. Ein komplettes Minimalbeispielmit allen erwähnten Funktionsaufrufen ist in Anhang A zu finden.

2.2.2. Device-Code / Kernel

Der Code, der auf den OpenCL-Devices ausgeführt wird, wird Kernel genannt.Ein Kernel ähnelt einer, zum Beispiel in C oder C++ üblichen Funktion, dievom Host aufgerufen wird. Ein Kernel stellt den Einstiegspunkt für die Aus-führung von Code auf dem Device dar.Ein kurzes Beispiel soll genauer erläutern, was einen Kernel von einer üblichenC-Routine unterscheidet und wie damit die extrem hohe Parallelität erreichtwird, für die OpenCL ausgelegt ist.

Bei großen Datensätzen ist es üblich, diese mittels Schleifen zu durchlaufen.Müssen zum Beispiel mehrdimensionale Daten bearbeitet werden, erfolgt diesüblicherweise mit verschachtelten Schleifen:

c© Maik Dankel 18


2. GPGPUfor (i=0; i


2. GPGPUdagegen ist eine spezielle Instanz des Kernels, die einen bestimmten Daten-satz bearbeitet. Angelehnt an das vorherige Beispiel entspricht ein Kernel derFunktion process(point[i][j][k]). Ein einzelnes Work-Item entspricht zum Bei-spiel process(point[1][2][3]).

Jedes Work-Item besitzt eine eindeutige ID, die Global-ID genannt wird. DieGlobal-ID der Kernelinstanz process(point[1][2][3]) ist das Array {1, 2, 3}.Die Anzahl an Elementen, aus der eine Global-ID besteht, steht in direktemBezug zur Dimension der Daten, die bearbeitet werden. Die Dimension wirdbeim Aufruf von clEnqueueNDRangeKernel mit dem Argument work-dim an-gegeben. Bei dem vorherigen Beispiel ist work-dim gleich drei gesetzt. Darausergibt sich folgender Kernel-Code:

int i = get_global_id(0);int j = get_global_id(1);int k = get_global_id(2);process(point[i][j][k]);

Der Aufruf von get_global_id(n) liefert die eindeutige ID, der n-ten Dimensi-on.

Das nächste Konzept, das in diesem Zusammenhang zu erwähnen ist, ist wasNVIDIA in seiner Architektur als CUDA Thread Blocks bezeichnet und inOpenCL Work-Group genannt wird. Eine Work-Group ist ein Zusammen-schluss mehrerer Work-Items, die gemeinsam auf den extrem schnellen lokalenSpeicher zugreifen (siehe dazu Unterabschnitt 2.2.3) und untereinander syn-chronisiert werden können (siehe Unterabschnitt 2.2.4).Zusätzlich zur Global-ID besitzt ein Work-Item eine Local-ID, die es eindeu-tig innerhalb einer Work-Group identifiziert und mit get_local_id(n) abgefragtwerden kann (siehe Abbildung 2.5). Work-Groups bieten den Vorteil, dass nureine gewisse Anzahl an Work-Items die selben Daten bearbeiten und nicht al-le. Sollen beispielsweise viele 8x8-Matrizen gleichzeitig berechnet werden, sowird beim Aufruf von clEnqueueNDRangeKernel für local_work_size {8, 8}angeben. Dies bewirkt, dass eine Work-Group aus 8*8 Work-Items besteht.

In Abbildung 2.6 ist dargestellt, wie die vorhergehenden Prinzipien an einemzweidimensionalen Beispiel aussehen. Die einzelnen Quadrate stellen die Work-Items dar, die zu Work-Groups zusammengefasst sind. Global Size(0) bezeich-net die Gesamtanzahl Work-Items in der ersten und Global Size(1) in der

c© Maik Dankel 20


2. GPGPU

Abbildung 2.6.: Veranschaulichung von Work-Groups, Work-Items undDimensionen

Quelle: [J. Tompson 2012]

zweiten Dimension. Äquivalent dazu steht Local Size(n) für die Elemente in-nerhalb der Work-Groups.

2.2.3. Speichermodell

Wie bereits in Abschnitt 2.1 erläutert, verfügt eine GPU über verschiedeneSpeicherarten, die unterschiedlich schnelle Zugriffszeiten haben. Dies spiegeltsich durch verschiedene Speichertypen in OpenCL wieder. Die Speicherhier-archie (siehe Abbildung 2.7) in OpenCL ähnelt der physikalischen Speicher-konfiguration von NVIDIA und ATI Hardware, auch wenn das Mapping nichteins zu eins ist, da die beiden Hersteller ihre Speicherhirarchie unterschiedlichdefinieren [J. Tompson 2012].

c© Maik Dankel 21


2. GPGPU

Abbildung 2.7.: OpenCL SpeichermodellQuelle: [KHRONOS 2008]

Wie in Abbildung 2.7 zu sehen, verfügt eine OpenCL Umgebung über folgendeSpeicherbereiche:

• Host MemoryAls Hostspeicher wird im Allgemeinen der Hauptspeicher, oder Massen-speicher wie Festplatten, des Hosts bezeichnet.

• Global MemoryJede Work-Group und damit auch alle Work-Items haben sowohl lesen-den, als auch schreibenden Zugriff auf den globalen Speicher (siehe L2-Cache aus Unterabschnitt 2.1.1). Die OpenCL Spezifikation trifft hierbeikeine Aussagen, ob die Zugriffe gecached werden. Dies hängt von derjeweiligen Hardware ab [M. Müller 2009].

• Constant MemoryEin Teil des globalen Speichers wird für konstante Werte genutzt, diezur Laufzeit unverändert bleiben. Dieser Bereich wird nur hostseitig be-schrieben.

c© Maik Dankel 22


2. GPGPU• Local Memory

Die einzelnen Work-Groups verfügen über lokalen Speicher, der inner-halb dieser erreichbar ist (siehe Shared Memory / L1-Cache aus Unter-abschnitt 2.1.2). Alle Work-Items einer Work-Group können auf diesenSpeicher zugreifen.

• Private MemoryJedes Work-Item verfügt zusätzlich über einen kleinen privaten Speicher,in dem private Variablen gespeichert werden.

Der Global Memory entspricht dem in Unterabschnitt 2.1.1 angegebenen L2-Cache der GPU. Der Local Memory und der Private Memory entsprechendem Shared Memory, bzw. L1-Cache eines SMX, wie in Unterabschnitt 2.1.2beschrieben.

Der Transfer von Daten muss immer explizit folgende Schritte durchlaufen:Zunächst werden Daten beim Kernelaufruf vom Hostspeicher in den globalenDevicespeicher geladen. Auf diesem kann direkt gearbeitet werden. Da dieserSpeicher allerdings langsamer ist als der lokale Speicher, sollten die Daten, diebearbeitet werden, zunächst in den lokalen Speicher kopiert werden. Nach derBearbeitung werden diese zurück in den globalen Speicher geschrieben und kön-nen vom Host ausgelesen werden. Die oben angegebene Reihenfolge entsprichtdemnach gleichzeitig der Speicherhierarchie, von Host Memory (langsam) bishin zu Private Memory (sehr schnell).Was beim Transfer von Daten ebenfalls zu beachten ist, ist dass oftmals diePCIe Schnittstelle der Flaschenhals bei der Performance einer OpenCL An-wendung ist. Aus diesem Grund sollten Kopiervorgänge zum und vom Devicemöglichst immer große Datenmengen beinhalten. Viele kleine Datentransfersbremsen die Anwendung aus.

2.2.4. Asynchrone Verarbeitung und Synchronisation

OpenCL ist auf hohe Parallelität ausgelegt. Beim Aufruf von clEnqueueN-DRangeKernel werden beliebig viele Kernelinstanzen gestartet, die parallelDaten bearbeiten. Das bewirkt bereits eine asynchrone Verarbeitung, die alsData Parallelism bezeichnet wird und die Kernkompetenz von OpenCL-Anwendungen darstellt. Viele Kernelinstanzen führen parallel den selben Codeauf verschiedenen Daten aus. Arbeiten viele tausend Threads auf den Daten

c© Maik Dankel 23


2. GPGPUwird es spätestens bei etwas komplexeren Anwendungen von Nöten sein, dassdiese sich untereinander synchronisieren, da es vorkommen kann, dass zur Be-arbeitung eines Datums eine vorherige Berechnung abgeschlossen sein muss.Um dies zu bewerkstelligen, sind in OpenCL sogenannten Barriers, also zuDeutsch Barrieren, implementiert, an denen die Work-Items ”warten” müssen,bis alle an dieser Stelle im Code angelangt sind.

Abbildung 2.8.: Synchronisierung von Work-ItemsQuelle: [J. Tompson 2012]

Solch eine Synchronisation ist (wie in Abbildung 2.8 dargestellt) nur innerhalbeiner Work-Group möglich. Das heißt Work-Items können nicht Work-Group-übergreifend synchronisiert werden.Ein Synchronisationspunkt wird in OpenCL mit dem in Listing 2.14 gezeigtenAufruf von barrier() erzeugt. Jedes Work-Item einer Work-Group muss dieseFunktion ausführen, bevor mit der weiteren Ausführung fortgefahren werdendarf. Zusätzlich werden mit diesem Aufruf Speicheroperationen synchronisiert.Über den Parameter flags kann angegeben werden, welche Speicherbereichedavon betroffen sind. Als möglicher Parameter existiert CLK_LOCAL_MEM_FENCE oder CLK_GLOBAL_MEM_FENCE. Diese geben ab ob Opera-tionen auf dem lokalen, oder globalen Speicher synchronisiert werden5.

1 void barrier ( cl_mem_fence_flags flags)

Listing 2.14: OpenCL Synchronisationspunkt barrier()

5die Parameter können durch ein logische ODER kombiniert werden, wodurch auf Speicher-operationen auf beiden Speicherarten reagiert wird

c© Maik Dankel 24


2. GPGPUZu beachten ist, dass diese Synchronisationspunkte im Code immer von allenWork-Items erreichbar sein müssen. Ein barrier() Aufruf sollte nicht in einemIf-Statement stehen, falls dies bewirkt, dass nicht alle Instanzen diesen Codeausführen. Dies könnte zu Fehlern oder Deadlocks führen, da die blockiertenThreads nicht weiterlaufen, bis alle den Synchronisationspunkt erreicht haben,einige diesen aber nicht erreichen können. Selbiges gilt für Schleifen.

Zusätzlich zum Data Parallelism unterstützt OpenCL das Konzept des TaskParallelism. Es können verschiedene Aufgaben (Kernel) gleichzeitig ausge-führt werden. Dies kann zum Beispiel gewünscht sein, wenn die selben Datenauf unterschiedliche Art und Weise verarbeitet werden müssen. Hierzu werdenmehrere verschiedene Kernel aufgerufen und gestartet.

Ebenfalls lässt sich dieses Prinzip einsetzen, um eine GPU möglichst effizientauszulasten. Wenn ein Kernel den Grafikprozessor nicht genug beanspruchtund ein Großteil der Zeit auf die Kopiervorgänge von und zu diesem gewartetwird, kann der selbe Kernel mehrfach mit unterschiedlichen Daten aufgerufenwerden. Dadurch wird höhere Auslastung und damit höhere Parallelität er-reicht.Soll dieses Prinzip angewandt werden, muss wie in Unterabschnitt 2.2.1 - ”Spei-cher allozieren” erwähnt, mit den von OpenCL zur Verfügung gestellten Eventsgearbeitet werden. Werden diese nicht verwendet, sind sämtliche clEnqueue*Aufrufe blockierend und es können nicht mehrere Kernel gleichzeitig gestartetund ausgeführt werden.

1 cl_int clWaitForEvents (2 cl_uint num_events,3 const cl_event ∗event_list)

Listing 2.15: clWaitForEvents() Signatur

Der Einsatz dieser Technik stellt den Programmierer vor ein erneutes Syn-chronisationsproblem. Hierfür sind die OpenCL-Events zuständig. Wird beimAufruf einer clEnqueue* Routine eine Referenz auf ein cl_event übergeben, soenthält dieses Informationen über die Ausführung dieser speziellen Kernelin-stanz. Mit dem Aufruf von clWaitForEvents (siehe Listing 2.16) kann hostseitigeine Synchronisation vorgenommen werden, da dieser Aufruf blockiert, bis dieAusführung der Kernelinstanz beendet ist. Das ist zum Beispiel dann wichtig,wenn der Host mit Werten rechnet, die vom Device vorverarbeitet werden.

c© Maik Dankel 25


2. GPGPUWird an solch einer Stelle keine Synchronisation vorgenommen, ist die Daten-integrität gefährdet.Mit clWaitForEvents wird auf ein oder mehrere spezielle Events gewartet. DasArgument num_events gibt die Anzahl an Events an, auf die gewartet wirdund *event_list eine Referenz auf die entsprechenden Events. Soll nicht aufspezielle Events gewartet werden, sondern auf die Beendigung aller bislang ineiner Command-Queue abgesetzten Befehle kann alternativ clFinish gerufenwerden.

1 cl_int clFinish (cl_command_queue command_queue)

Listing 2.16: clFinish() Signatur

Hier wird die entsprechende Command-Queue übergeben und die weitere Aus-führung des Host-Codes solange blockiert, bis die Warteschlange deviceseitigabgearbeitet ist.

c© Maik Dankel 26


3. Kalman-Filter

3. Kalman-Filter

Der Kalman-Filter ist eine nach seinem Entdecker Rudolf Emil Kálmán be-nannte Sammlung von mathematischen Gleichungen. Dieser veröffentlichteerstmals 1960 in seinem berühmten Paper ”A New Approach to Linear Filte-ring and Prediction Problems” [R. E. Kalman 1960] eine rekursive Lösung zurlinearen Filterung von diskreten Daten. Vereinfacht gesagt erlaubt der Kalman-Filter bei Vorliegen von fehlerbehafteten Beobachtungen (zum Beispiel Mes-sungen) Rückschlüsse auf das System zu ziehen und damit Fehler zu minimie-ren. Mathematisch gesehen entspricht der Kalman-Filter einem Bayes’schenMinimum-Varianz-Schätzer für lineare stochastische Systeme.

Der Kalman-Filter ist ein Satz von mathematischen Gleichungen6, deren Be-sonderheit ist, dass sie besonders effizient in Echtzeitsystemen einsetzbar sind,um den Zustand eines Prozesses zu ”schätzen”. Aus diesem Grund findet er seitjeher rege Beachtung in der Wissenschaft und Wirtschaft. Dazu zählen unteranderem Anwendungen zur autonomen Navigation, Auswertung von Radar-signalen oder die Positionsverfolgung sich bewegender Objekte (Tracking)[G.Welch 2006].

In dieser Arbeit wird ein Kalman-Filter eingesetzt, um die Genauigkeit derMessungen aus dem Atlas-Teilchendetektor zu erhöhen und damit dem tat-sächlichen Punkt, den dieses Teilchen passiert, näher zu kommen. Das Prinzipdieses Algorithmus wird nachfolgend am Beispiel einer Spurrekonstruktion ge-nauer erläutert.

3.1. Definition

Die rekursive Form des Kalman-Filters kommt dadurch zustande, dass er einenkontinuierlichen Prozess darstellt, der sich laufend wiederholt, solange neue

6Die genaue Herleitung der Formeln ist nicht Teil dieser Arbeit und kann z.B. in [R. E.Kalman 1960] oder [G. Welch 2006] eingesehen werden

c© Maik Dankel 27


3. Kalman-FilterMesswerte vorhanden sind. Dieser Algorithmus hat einen erheblichen Vorteilgegenüber vergleichbaren Ansätzen, da er zur Berechnung eines Schrittes le-diglich die vorherige Berechnung benötigt, nicht aber alle vorangegangenenDaten (wie es zum Beispiel beim Wiener-Filter der Fall ist)[G. Welch 2006].Daher ist der Kalman-Filter sehr effizient, sowohl was Geschwindigkeit, alsauch Speicherbedarf angeht, ohne dass Informationen verloren gehen, da einberechneter Wert immer die Information aller vorhergehender Werte beinhal-tet.

Abbildung 3.1.: Der rekursive Ablauf des Kalman-Filters. Der time updateSchritt bildet den aktuellen Zustand auf einen geschätzten zu-künftigen Wert ab. Der measurement update Schritt korrigiertden geschätzten Wert mittels einer aktuellen Messung.

Quelle: [G. Welch 2006]

Abbildung 3.2 veranschaulicht schematisch, wie ein Teilchen verschiedene De-tektorlagen durchdringt und dabei Messpunkte erzeugt. Allerdings entsprechendiese Messungen nicht immer genau der Realität. Sie werden durch verschie-dene Faktoren beeinflusst, wie beispielsweise der Auflösung der jeweiligen De-tektorlage, Magnetfelder innerhalb des Detektors und anderer Störungen. Ausdiesem Grund wird der Kalman-Filter eingesetzt, um mittels Minimierung derFehlerkovarianz die fehlerbehaftete Messung zu korrigieren und somit näheran die tatsächlichen Durchschlagpunkte eines Teilchens zu gelangen.Der Kalman-Filter kann, wie in Abbildung 3.1 dargestellt, in zwei Schritte un-terteilt werden: der Time Update Schritt und derMeasurement Uptdate Schritt.Ersterer stellt eine Vorhersage (Schätzung) eines Zustandes dar, welcher imzweiten Schritt mittels einer konkreten Messung korrigiert wird.Im Time Update Schritt erfolgt zunächst eine Projektion von Lage k-1 zu La-ge k. Hierzu wird der vorherige Status-Vektor pk−1|k−1 mit der zugehörigen

c© Maik Dankel 28


3. Kalman-Filter

Abbildung 3.2.: Schematische Darstellung von Teilchenspur und verfälschterMessung innerhalb eines Detektors

Jacobi-Matrix Fk multipliziert und somit eine a-priori-Schätzung des nächstenSchritts berechnet. Daraus ergibt sich nachfolgende Gleichung 3.1

pk|k−1 = Fkpk−1|k−1 (3.1)

Neben dem Status (welcher im Beispiel Spurrekonstruktion einen Punkt imDetektor darstellt) wird auch für die zugehörige Fehlerkovarianzmatrix Ck|k−1eine a-priori-Schätzung berechnet. Diese besteht aus der Summe von zwei Ter-men, die in Gleichung 3.2 angegeben sind. Der erste beschreibt die lineare Feh-lerfortpflanzung, der vorherigen Fehlerkovarianzen, der zweite einen Rausch-term.

Ck|k−1 = FkCk−1|k−1F Tk + PkQkPTk (3.2)

Nach Berechnung von Gleichung 3.1 und Gleichung 3.2 ist die Time UpdatePhase abgeschlossen und es folgt die sogenannte Measurement Update Phase.Hier werden die zuvor geschätzten Ergebnisse von pk|k−1 und Ck|k−1 unter Be-rücksichtigung eines neuen Messwertes korrigiert und aktualisiert und somitein a-posteriori-Wert berechnet.Dazu wird der sogenannte Kalman-Gain Kk (siehe Gleichung 3.3) berechnet,der angibt, wie stark diese Korrektur in den neuen Wert einfließt. Die Ma-trix Hk ist hierbei und im Folgenden eine Transformationsmatrix, die dafür

c© Maik Dankel 29


3. Kalman-Filtersorgt, dass die einzelnen Parameter die korrekte Dimensionalität haben, Vkbeschreibt einen Rauschterm.

Kk = Ck|k−1HTk (Vk + HkCk|k−1HTk )−1 (3.3)

Gleichung 3.4 stellt die Korrektur des zuvor geschätzten Wertes von pk|k−1unter Berücksichtigung eines aktuellen Messwertes mk dar. Es wird folglicheine a posteriori Abschätzung des Statusvektors berechnet.

pk|k = pk|k−1 + Kk(mk −Hkpk|k−1) (3.4)

Abschließend wird mittels Gleichung 3.5 die korrigierte FehlerkovarianzmatrixCk|k berechnet.

Ck|k = (I −KkHk)Ck|k−1 (3.5)

Gleichung 3.1 und Gleichung 3.2 stellen den Time Update Schritt dar, Glei-chung 3.3 bis Gleichung 3.5 entsprechen dem Measurement Update Teil desKalman-Filters7.Nach Beendigung beider Schritte werden die aktuellen a posteriori Schätzun-gen von pk|k und Ck|k als Input für den nächsten Time Update Schritt ver-wendet und der in Abbildung 3.1 dargestellte Zyklus startet von vorn. Abbil-dung 3.3 veranschaulicht schematisch die durch den Filter optimierte Spur, dienäher an den tatsächlichen Werten liegt.

3.2. Startwerte

Eine gewisse Herausforderung stellt das Finden von geeigneten Startwertenfür den Filter dar. Im ersten Schritt existiert kein ”Vorgänger” mit Index k-1auf den zugegriffen werden kann. Werden schlechte Startwerte gewählt, kön-nen diverse Probleme auftreten. Die lineare Approximation des Spurmodellskönnte falsch sein, falls der Startwert zu weit vom tatsächlichen Wert ent-fernt liegt. Außerdem könnte die Vorhersage gänzlich fehlschlagen, weil derberechnete Pfad die nächste Detektorlage nicht schneidet. Ein Ansatz zur Lö-sung dieser Problematik ist eine Referenzspur zu verwenden, die zuvor durch

7Die angegebenen Gleichungen variieren in verschiedenen Publikationen leicht. Die hierverwendete Notation stammt aus [R. Frühwirth 2000], S. 246 f.

c© Maik Dankel 30


3. Kalman-Filter

Abbildung 3.3.: Schematische Darstellung der durch den Kalman-Filter korri-gierten Spur

Mustererkennungsverfahren ermittelt wird. Es werden zunächst für die Refe-renzspur sämtliche Schnittpunkte mit den Detektorlagen berechnet und alsReferenzwerte pk,r verwendet. Der Filterprozess findet nicht mehr auf den rei-nen Messwerten statt, sondern auf der Differenz zu dieser zur Referenzspur.Daraus ergibt sich, dass nicht mk, sondern ∆mk = mk −Hkpk,r verwendetwird. Dieses Prinzip verbessert die Qualität des Filter-Outputs und erleich-tert das Finden geeigneter Startwerte ([R. Frühwirth 2000], S. 252). Da eineDifferenz hier im Optimalfall null beträgt, kann ohne weiteres Vorwissen diesals Startwert genutzt werden. Wird keine Referenzspur genutzt und null alsStartwert gewählt, kann die Abweichung zum eigentlichen Messpunkt beliebiggroß werden. Daher wird empfohlen, wann immer möglich eine Referenzspur zuverwenden. Dieses Prinzip auf der Differenz zu einer Referenzspur zu arbeitenwird ebenfalls in dieser Arbeit verwendet (siehe Unterabschnitt 4.3.1).

3.3. Smoothing

Der Kalman-Filter liefert sehr gute Ergebnisse. Die in Abschnitt 3.1 beschriebe-nen Arbeitsweise des Filters, bei der ein Wert immer nur von seinem Vorgängerabhängt, hat jedoch einen nicht zu verachtenden Nachteil. Jeder neu berech-nete Wert enthält immer sämtliche Informationen der bis dahin errechnetenWerte, da diese aufeinander aufbauen. Das bedeutet folglich, dass, je weiter

c© Maik Dankel 31


3. Kalman-Filterder Prozess fortgeschritten ist, der Informationsgehalt und damit die Wahr-scheinlichkeit genauere Werte zu erhalten steigt. Im Umkehrschluss bedeutetdies aber, dass am Anfang der Reihe sehr wenig Information vorhanden istund die Wahrscheinlichkeit, dass die errechneten Werte ungenau sind, relativhoch ist.Um diesem Phänomen entgegenzuwirken, kann eine Technik namens Smoo-thing (zu deutsch: Glättung) eingesetzt werden. Ziel ist hierbei diesem Un-gleichgewicht des Informationsgehaltes entgegenzuwirken. Dazu kann der Kal-man-Filter noch einmal auf den selben Daten, jedoch in entgegensetzter Rich-tung ausgeführt werden. Das bewirkt, dass die Informationsgehalte der einzel-nen Zwischenergebnisse in entgegengesetzter Richtung zunehmen. Werden dieentsprechenden Werte pf und pb (das hochgestellte ’f’ steht für forward, alsoden Hinweg, das ’b’ entsprechend für backward, also den Rückweg) zusammen-geführt und gewichtet, wird dadurch ein konsistent hoher Informationsgehaltüber die gesamte Messreihe erreicht.Gleichung 3.6 gibt an, wie die Gewichtung der einzelnen Messwerte berech-net wird. Hierzu werden die Fehlerkovarianzmatrizen aus Hin- und Rückwegverwendet.

K ′k = Cfk|k ∗ (C

fk|k + C

bk|k)−1 (3.6)

Gleichung 3.7 stellt die Berechnung des geglätteten Wertes dar. Zu beach-ten ist, dass aus dem Rückwärts-Fit der a-priori-Wert pbk|k−1 und aus demVorwärts-Fit der a-posteriori-Wert pfk|k verwendet wird, da sonst die jeweiligeMessung im Smoothing effektiv doppelt gewichtet würde.

p′k = pfk|k + K

′k ∗ (p

bk|k−1 − p

fk|k) (3.7)

Letztendlich kann zu dem neuen, geglätteten state-Vektor eine zugehörige Feh-lerkovarianzmatrix C ′k berechnet werden. Die Berechnung ist in Gleichung 3.8dargestellt.

C ′k = (I −K′k)C

fk|k (3.8)

c© Maik Dankel 32


4. Implementierung

4. Implementierung

Wie eingangs erwähnt, ist Ziel dieser Arbeit die Implementierung eines, wie inKapitel 3 beschriebenen, Kalman-Filters zur Rekonstruktion der Spur, die einTeilchen zurücklegt, wenn es sich durch die einzelnen Detektorlagen bewegt.Da extrem große Datenmengen entstehen, wenn zwei Teilchenströme im Be-reich eines Detektors aufeinandergelenkt und zur Kollision gebracht werden,ist das Ziel die Auswertung dieser Daten möglichst mit einem hohen Gradan Parallelität vorzunehmen. Um dies zu erreichen wird das in Abschnitt 2.2erläuterte OpenCL-Framwork genutzt. Nachfolgend wird das Zusammenspielund speziell die Umsetzung der zuvor erläuterten Technologien und Algorith-men genauer beschrieben.

4.1. Data-IO und Datenstrukturen

Die auszuwertenden Testdaten liegen in diesem Fall in Form eines ROOT-Filesvor. ROOT ist eine 1995 vom CERN entwickelte, objektorientierte Softwa-re, die insbesondere mit großen Datenmengen umgehen und diese analysierenkann. Es ermöglicht einem (sowohl Endanwender, als auch Entwickler, der die-ses Framwork nutzt) Daten effizient und ohne diese verändern zu müssen, zukombinieren und zu visualisieren. So stehen einem Funktionalitäten wie Hi-stogrammierung, Curve Fitting, Auswertung von Funktionen, sowie weitereGraphik- und Visulisierungsklassen zur Verfügung.Der hier verwendete Datensatz liegt in Form einer Baumstruktur (TTree) vor,die in einer Klasse EventReader ausgelesen und in Events (Kollisionen vonTeilchen) und zugehörige Tracks (Spuren der Teilchen) gespeichert werden.8

Listing 4.1 zeigt solch ein Event, welches der EventReader ausliest. Diesesbesteht wiederum aus mehreren Tracks, weshalb es als Vektor von Tracks mo-delliert ist.

8Anm.: Die Implementierung der Klasse EventReader war Teil eines Masterprojektes undnicht Bestandteil dieser Arbeit

c© Maik Dankel 33


4. Implementierung

1 typedef std :: vector KF_Event_t;

Listing 4.1: Definition eines Events

Die Track_t-Struktur wird in diesem Fall so modelliert, dass sie aus mehrerenInformationen besteht. Zum einen dem eigentlichen Track, der wiederum ausmehreren Hits (also Messpunkten im Detektor) besteht und hier mit trackbetitelt ist. Zum anderen aus dem Startpunkt der Referenzspur (info) unddem tatsächlichen Startpunkt der Simulation (truthTrackInfo).

1 typedef std :: vector TrackData_t;2

3 struct TrackStruct {4 TrackData_t track;5 TrackInfo_t info;6 TrackInfo_t truthTrackInfo;7 };8 typedef TrackStruct Track_t;

Listing 4.2: Definition eines Tracks

Die Daten, die ein Track beinhaltet (TrackData_t) sind, wie in Listing 4.2gezeigt, als Vektor von Trackhits (TrackHit_t) modelliert. Der Aufbau einesTrackhits ist in Listing 4.3 abgebildet.

1 struct TrackHitStruct {2 scalar_t normal[ORDER];3 scalar_t ref [ORDER];4 scalar_t err_locX;5 scalar_t err_locY;6 scalar_t cov_locXY;7 scalar_t jacobi [ORDER ∗ ORDER];8 scalar_t jacobiInverse [ORDER ∗ ORDER];9 char is2Dim;

10 int detType;11 int bec;12 };13 typedef struct TrackHitStruct TrackHit_t;

Listing 4.3: Definition eines Tracks

c© Maik Dankel 34


4. ImplementierungWie in Listing 4.3 zu sehen, beinhaltet ein einzelner Hit sehr viele Informa-tionen. So existiert zu jedem Hit ein Array ref in der Dimension ”ORDER”(ORDER ist ein Projektweites #define und ist in diesem speziellen Fall auf 5festgelegt), das einen Vektor darstellt, der den Punkt der Refernzspur im De-tektor repräsentiert. Der Datentyp dieser Werte ist scalar_t, welcher über ein#define auf float oder double (oder sonstige benötigte Datentypen) festgelegtwerden kann. Der Vektor sieht folgendermaßen aus:

~ref =

refLocX

refLocY

refPhi

refTheta

refQoverP

refLocX und refLocY geben die Position des Hits in X/Y-Koordinaten an, ref-Phi und refTheta beschreiben Winkel der Flugbahn und refQoverP beschreibtdie Krümmung der Flugbahn. Entsprechend dazu existiert ein Vektor der si-mulierten Hits, bei dem allerdings die Werte Phi, Theta und QoverP gleich nullsind, da ein Messwert nur über einen Punkt im Koordinatensystem verfügt,nicht aber über die Eigenschaften eines sich bewegenden Teilchens.

~normal =

locX

locY

000

Des Weiteren existiert zu einem Hit eine Matrix, die die Propagation von De-tektorlage i-1 nach Lage i beschreiben. Das ist die sogenannte Jacobi-MatrixFk, die in Gleichung 3.1 und Gleichung 3.2 auftaucht und hier mit jaco-bi[ORDER*ORDER] bezeichnet ist. Der Einfachheit halber wird in dieser Da-tenstruktur zusätzlich die Inverse dieser Matrix unter jacobiInverse abgelegt,da diese beim Rückwärtsfit benötigt wird.Außerdem existiert zu jedem Hit eine Information über den Detektortyp (det-Type), der angibt, in welchem Bereich des Detektors sich der Hit befindet.Zusätzlich ist ein Wert is2Dim vorhanden, welcher ein Flag darstellt, das an-gibt, ob die Messung eindimensional, oder zweidimensional berücksichtigt wer-den soll. Dies stellt einen Spezialfall dar, welcher in Abschnitt 4.2 genauer er-

c© Maik Dankel 35


4. Implementierungläutert wird. Die Werte err_locX, err_locY und cov_locXY werden für dieMatrix Vk aus Gleichung 3.3 benötigt.

4.2. Host-Code

Abbildung 4.1.: Arbeitsweise der Implementation

Der Ablauf der implementierten Software kann, wie in Abbildung 4.1 darge-stellt, grob in sechs Stufen eingeteilt werden. Vier davon sind Aufgaben, dieauf dem Host-System ausgeführt werden (blau hinterlegt) und zwei Schritte er-folgen auf dem Device (rot hinterlegt). Nachfolgend wird die hostseitige Imple-mentierung erläutert und in Abschnitt 4.3 der entsprechende Device-Code.

Im ersten Schritt werden, unter Verwendung des EventReaders, Daten aus demROOT-File ausgelesen und in den in Listing 4.1 und Listing 4.2 angegeben Da-tenstrukturen gespeichert.Anschließend werden die für die Ausführung des OpenCL-Codes erforderlichenElemente, wie program und kernel, erstellt und initialisiert. Im Anschluss star-tet der Hauptteil der Software in Form einer Schleife, die alle Events und alledarin enthaltenen Tracks bearbeitet.Wie im Pseudocode in Listing 4.4 zu sehen ist, erfolgt die Bearbeitung derDaten auf Seiten des Hosts seriell. Es werden sowohl die Events, als auch diedarin enthaltenen Tracks nacheinander abgearbeitet. Es folgt zunächst eine In-titalisierungsphase, bei der die zur Verfügung stehenden Daten vorverarbeitetwerden und Speicherplatz reserviert wird.Abbildung 4.2 stellt schematisch dar, wie die Aufteilung des Speichers erfolgt,in dem die zu bearbeitenden Daten gespeichert werden. Es werden für dieWerte, die vom Kalman-Filter berechnet werden (fits[]), (ORDER*hitcount)*3

c© Maik Dankel 36


4. ImplementierungElemente, also dreimal die Anzahl der Hits pro Track, reserviert. Grund da-für ist, dass in diesem Array zunächst die Werte aus dem Vorwärts-Fit, danndie Werte des Rückwärts-Fit und zuletzt die Ergebnisse aus dem Smoothingabgelegt werden.

1 for ( i = 0; i < #Events; i++){2 for (j = 0; j < #Tracks; j++) {3 init ()4 kernelFiltering ()5 MatrixInvert()6 kernelSmoothing()7 handleResults()8 }9 }

Listing 4.4: Pseudocode Filtering-Prozess

Der Speicherbereich für die Kovarianzmatrizen C_kk ist so angelegt, dass diejeweiligen Matrizen aus dem Vorwärts- und dem Rückwärtsfit abgespeichertwerden können. Da es sich hierbei um 5x5 Matrizen und nicht Vektoren han-delt, ist hierbei der Faktor ORDER2 angegeben.Diese Anordnung wird bewusst so gewählt, da dies die Daten sind, die auf dem

Abbildung 4.2.: Speicheraufteilung

OpenCL-Device berechnet werden. Das bedeutet, dass diese Daten vom undzum Device transferiert werden müssen. Da, wie in Unterabschnitt 2.2.3 be-reits erwähnt, viele kleine Datentransfers nach Möglichkeit vermieden werdensollen, wird so die Möglichkeit geschaffen, einen möglichst großen Datenbereichkomplett zu kopieren.

c© Maik Dankel 37


4. Implementierung

1 int hitCount = reader.getTrack(i, j) .track. size () ;2 TrackHit_t hits[hitCount];3 for ( int k = 0; k < hitCount; k++) {4 hits [k] = reader.getTrack(i, j) .track.at(k);5 if (abs(hits [k ]. normal[1] + 99999.) > 1.E−4)6 hits [k ]. is2Dim[k] = true;7 else8 hits [k ]. is2Dim[k] = false ;9 }

Listing 4.5: Überprüfung der Dimensionalität

Eine weitere Besonderheit, die es zu beachten gilt, ist die Behandlung des inAbschnitt 4.1 erwähnten Spezialfalls. Im Allgemeinen werden die Messungenaus dem Detektor zweidimensional berücksichtigt. Das heißt, ein Messwert be-steht aus einer X- und einer Y-Koordinate. Auf Grund der Detektorarchitekturist dies allerdings nicht immer der Fall, da der Silizium-Streifendetektor (SCT,Semiconductor Tracker) nur eine X-Koordinate besitzt. In diesem Fall wirdder Default-Wert ”-99999” für die Y-Koordinate gesetzt. Dies muss im Codeüberprüft und abgefangen werden, da es sonst zu erheblich verfälschten Fit-Resultaten führt.In Listing 4.5 ist der entsprechende Code-Abschnitt dargestellt. Es wird zu-nächst die Anzahl der Hits pro Track mit Hilfe des EventReaders bestimmtund in einem Integer-Wert hitCount gespeichert. Danach wird in einer Schleifedas Array hits mit den Hits gefüllt und gleichzeitig die zuvor erwähnte Über-prüfung durchgeführt und das Flag is2Dim entsprechend gesetzt.

Nach Durchführung dieser Schritte, wird Speicherplatz auf dem OpenCL -Device reserviert. Dies geschieht unter Verwendung des in Abschnitt 2.2.1 be-schriebenen Aufrufs von clCreateBuffer. Zum Device transferiert wird die zuvorintialisierte Datenstruktur hits, sowie der Wert hitCount. Für die Arrays fitsund C_kk wird ein Speicherbereich mit Lese- und Schreibrechten reserviertund keine Daten kopiert, da diese erst auf dem Device berechnet und der Spei-cher gefüllt wird.

Damit ist die Initialisierung abgeschlossen und der Kernel doFiltering wirdmittels Aufruf von clEnqueueNDRangeKernel gestartet und auf dem Deviceausgeführt (siehe dazu Abschnitt 4.3).

c© Maik Dankel 38


4. Implementierung

1 scalar_t ∗invC = C_kk;2 invC += ORDER ∗ ORDER ∗ hitCount;3 for ( int i = 0; i < hitCount; i++) {4 reader.invertMatrix(invC, invC, ORDER, ORDER);5 invC += ORDER ∗ ORDER;6 }

Listing 4.6: Invertierung der Kovarianzmatrizen

Nach Ausführung des Kernels wird zunächst, wie in Listing 4.6 gezeigt, derSpeicherbereich, in dem die Kovarianzmatrizen gespeichert sind (C_kk), vomDevice zum Host kopiert. Hostseitig wird anschließend die in Gleichung 3.6angebene Invertierung (Cfk|k + Cbk|k)−1 vorgenommen9.Listing 4.6 zeigt den Codeabschnitt, in dem diese Invertierung durchgeführtwird. Es wird zunächst ein Zeiger auf den hinteren Teil des Speicherbereichsgesetzt (Vgl. Abbildung 4.2). Hierbei ist zu erwähnen, dass deviceseitig indiesen Bereich bereits Cfk|k + Cbk|k, also die ohnehin benötigte Summe derMatritzen aus Hin- und Rückweg, gespeichert wird. Dies ermöglicht direkt aufdiesem Speicherbereich die Invertierung durchzuführen. Die Invertierung selbstgeschieht unter Verwendung der GSL-Bibliothek und basiert auf der sogenann-ten Singular Value Decomposition (SVD) und berechnet die Pseudoinverse derangegebenen Matrix.

Nach der Invertierung wird der Speicherbereich, der die Inverse der Summeder Kovarianzmatrizen enthält, zurück auf das Device kopiert und der KerneldoSmoothing gestartet. Hier wird das Smoothing durchgeführt und anschlie-ßend die Endergebnisse zurück zum Host kopiert. Damit ist die hostseitigeImplementierung des Kalman-Filters abgeschlossen und die Ergebnisse kön-nen weiter verarbeitet werden.

4.3. Device-Code

Die Implementierung des Kalman-Filters findet, bis auf den zuvor beschriebe-nen Teil der Invertierung, ausschließlich deviceseitig statt. Dazu werden zweiOpenCL-Kernel entwickelt, die den in Abschnitt 3.1 erläuterten Kalman-Filter,

9Tests haben ergeben, dass die hostseitige Invertierung trotzdem erhöhtem Datentransferzwischen Host und Device zeiteffizienter ist

c© Maik Dankel 39


4. Implementierungsowie das in Abschnitt 3.3 beschriebene Smoothing auf einem OpenCL-Deviceausführen. Zusätzlich existieren einige Hilfsfunktionen, die die Les- und Wart-barkeit des Codes erhöhen und teilweise von beiden Kernels benötigt werden.

4.3.1. Kalman-Filter

Die Software bearbeitet hostseitig seriell Track für Track aus den einzelnenEvents. Es werden die Daten von jedem Track aufbereitet und auf das Devicekopiert. An dieser Stelle kommt die Fähigkeit von OpenCL-Anwendungen zutragen, Daten parallel verarbeiten zu können. Wie zuvor beschrieben, beste-hen die zu verarbeitenden Daten zumeist aus 5x5 Matrizen, welche als Arraysmit je 25 Elementen modelliert sind. Da der Kalman-Filter zu großen Teilenmit Matrixoperationen arbeitet, werden 25 Threads pro Spurrekonstruktionverwendet. Diese Threads sind in zwei Dimensionen aufgeteilt, sodass entspre-chend der Datenbeschaffenheit 5x5 Threads (die eine Work-Group bilden) zurVerfügung stehen. Jeder Thread berechnet eine Position in der Ergebnismatrix.Hier wird die hohe Parallelisierbarkeit von Matrixmultiplikationen ausgenutzt,um den Code zu beschleunigen.

1 __kernel void doFiltering (...) {2 init ()3 for ( i = 0; i < #Hits; i++){4 timeUpdate()5 measurementUpdate()6 }7 reinit ()8 for ( i = #Hits−1; i >=0; i−−){9 timeUpdate()

10 measurementUpdate()11 }12 }

Listing 4.7: Pseudo-Code Filtering-Kernel

Der Ablauf des Filterprozesses auf dem Device kann im Pseudo-Code wie inListing 4.7 dargestellt werden: Er besteht im Wesentlichen aus zwei Schleifen,die nacheinander abgearbeitet werden. Die Schleifen selbst stellen den Filter-vorgang dar. Die erste der beiden berechnet den Vorwärts-Fit und die zweiteden Rückwärts-Fit. Die Besonderheit hierbei ist, dass der Code innerhalb der

c© Maik Dankel 40


4. ImplementierungSchleifen gleichzeitig von den zuvor erwähnten 25 Threads ausgeführt wird.Die Schleifen selbst müssen allerdings von den jeweiligen Threads auf Grundder rekursiven Natur des Kalman-Filters seriell abgearbeitet werden und kön-nen nicht weiter parallelisiert werden.Wie in Unterabschnitt 2.2.2 erläutert können diese Threads über den Aufrufvon get_local_id() innerhalb einer Work-Group eindeutig identifiziert werden.Deshalb werden zunächst zwei Variablen column und row im privaten Spei-cher des Threads angelegt und mit der jeweiligen local-ID initialisiert. Überdiese wird anschließend die jeweilige Spalte und Zeile, in der der Thread arbei-tet, identifiziert. Bei dem hier gewählten Modell einer 5x5 Work-Group habencolumn und row also immer einen Wert zwischen null und vier.

1 __kernel void doFiltering(__global const TrackHit_t ∗hits,2 __global const int ∗hitCount, __global trackHitData ∗fits,3 __global scalar_t ∗C_kk) {4 int column = get_local_id(1);5 int row = get_local_id(0);6 int currentHitIndex;7 int privHitCount = hitCount[0];8

9 __local scalar_t jacobi[ORDER ∗ ORDER];10 __local scalar_t normalData[ORDER];11 __local scalar_t refData[ORDER];12 ...

Listing 4.8: Funktionskopf des doFiltering Kernels

Nachfolgend werden zwei weitere private Variablen currentHitIndex und priv-HitCount angelegt, wobei Erstere den Schleifenzähler für die Hauptschleifendarstellt und Letztere die übergebene Anzahl der Hits dieses Tracks speichert.Anschließend werden diverse Variablen im lokalen Shared-Memory erstellt, diefür die Berechnung des Kalman-Filters benötigt werden. Zeile 4 bis 7 in Lis-ting 4.8 zeigt die Deklaration der vier privaten Variablen, welche sich dadurchkennzeichnen kein Prefix zu besitzen. Im Gegensatz dazu weisen die Variablenim lokalen Speicher das Prefix __local vor. Da dieser jedoch nur immer voneinem Work-Item erreichbar ist und nur eine relativ geringe Größe besitzt,werden die Variablen, die innerhalb der gesamten Work-Group verfügbar seinmüssen, im lokalen Speicher angelegt. Die vier privaten Variablen jedoch wer-den sehr oft (meist zur Indizierung) verwendet, allerdings von jedem Thread

c© Maik Dankel 41


4. Implementierungunabhängig und können daher im privaten Speicherbereich angelegt werden,was die Zugriffszeit verringert.

Nach der Deklaration folgt die Initialisierung der, für den Kalman-Filter wich-tigen, Variablen. Es werden für den Filter folgende Startwerte festgelegt:

Ck−1|k−1 =

250 0 0 0 00 250 0 0 00 0 0.25 0 00 0 0 0.25 00 0 0 0 1E − 6

pk−1|k−1 =

00000

Da der Fit, wie in Abschnitt 3.2 erläutert, auf der Differenz zu einer Refe-renzspur erfolgt, wird als Startwert für den State-Vektor ~0 angenommen, wasbedeutet, dass keine Differenz zur Referenzspur vorhanden ist. Da dies jedochunwahrscheinlich ist, werden für die Kovarianzmatrix große Werte angegeben,was eine große Ungenauigkeit beschreibt.

Nach Initialisierung der Startwerte und Deklaration aller nötigen Variablenwird die iterative Ausführung des Vorwärts-Fits gestartet. Zu Beginn werdendie beim Aufruf des Kernels übergebenen Informationen des aktuellen Hits indie entsprechenden __local-Variablen geschrieben, um nicht auf dem langsa-meren globalen Speicher rechnen zu müssen.

1 for (currentHitIndex = 0; currentHitIndex < privHitCount; currentHitIndex++) {2 jacobi [row ∗ ORDER + column] = hits[currentHitIndex].jacobi[row ∗ ORDER+

column];3 normalData[row] = hits[currentHitIndex].normal[row];4 refData[row] = hits[currentHitIndex]. ref [row];5 err_locX = hits[currentHitIndex].err_locX;6 err_locY = hits[currentHitIndex].err_locY;7 cov_locXY = hits[currentHitIndex].cov_locXY;8

9 barrier (CLK_LOCAL_MEM_FENCE);10 ...11 }

Listing 4.9: Start des Vorwärts-Fits

c© Maik Dankel 42


4. ImplementierungWas in Listing 4.9 auffällt ist die Art und Weise der Indizierung, beispielsweiseder Jacobi-Matrix. Wie bereits erwähnt, ist die Position im Array, an demdas jeweilige Work-Item die Daten bearbeitet, eindeutig über column und rowbeschrieben. Dadurch wird mit der Anweisung

jacobi[row * ORDER + column] = ...

aus Zeile 2 in Listing 4.9 genau eine bestimmte Stelle des Arrays gefüllt, undzwar gleichzeitig für alle 25 Stellen. Ähnlich verhält es sich mit normalDataund refData. Dies sind jedoch Arrays der Größe 5 und werden daher nur mitdem Index row angesprochen. Das hat zur Folge, dass jedes Work-Item, dasdie selbe row-ID hat den Wert mehrmals an die selbe Stelle schreibt. Das stelltzwar eine Redundanz dar, eine Abfrage, ob der Wert bereits gesetzt ist, würdeallerdings erhebliche Performance-Einbußen zur Folge haben. Daher wird dieseRedundanz in Kauf genommen. Um sicherzustellen, dass alle Work-Items denSchreibvorgang auf dem lokalen Shared Memory abgeschlossen haben, wird ab-schließend eine Synchronisation mittels des barrier()-Aufrufs vorgenommen.

Anschließend startet die Time-Update Phase des Kalman-Filters. Hierzu wirdeine Hilfsfunktion mit Namen doPrediction() implementiert, um den Codeübersichtlicher und wartbarer zu halten. Diese Hilfsfunktion ist ebenfalls alsKernel implementiert, jedoch werden ihm Argumente aus dem lokalen SharedMemory übergeben und er kann daher nur innerhalb des eigentlichen Ker-nelcodes aufgerufen werden. Eine Besonderheit ist, dass diese internen Kernelimmer als Inline-Funktionen behandelt werden, was Kontextwechsel vermeidetund daher keinen Performance-Verlust zur Folge hat.

1 __kernel void doPrediction(__local scalar_t ∗jacobi, __local scalar_t ∗p_k_1k_1,2 __local scalar_t ∗p_kk_1, __local scalar_t ∗C_k_1k_1,3 __local scalar_t ∗C_kk_1) {4 //Predict state5 multiply5x1c_Av(jacobi, p_k_1k_1, p_kk_1);6 //Predict Error Covariance7 predict_C_kk_1(jacobi, C_k_1k_1, C_kk_1);8 }

Listing 4.10: Time-Update Hilfsfunktion

Der erste Schritt der Time-Update Phase entspricht nach Gleichung 3.1 ei-ner einfachen Matrix-Vektor-Multiplikation. Um die Multiplikation einer 5x5

c© Maik Dankel 43


4. ImplementierungTabelle 4.1.: Hilfsfunktionen Matrix-/Vektor-Multiplikation

Funktionsname Funktion Ergebnismultiply5x5 C = AB 5 x 5multiply5x5A_AB A = AB 5 x 5multiply5x5B_AB B = AB 5 x 5multiply5x1c_Av ~c = A~v 5 x 1multiply5x1v_Av ~v = A~v 5 x 1

Matrix mit einem 5x1 Vektor durchzuführen wird eine weitere Hilfsfunktioneingeführt, die nach dem Prinzip ~c = A~v arbeitet.

1 __kernel void multiply5x1c_Av(__local scalar_t ∗A, __local scalar_t ∗v,2 __local scalar_t ∗c) {3 int row = get_local_id(0);4 scalar_t result = 0;5 for ( int i = 0; i < ORDER; i++) {6 result += A[row ∗ ORDER + i] ∗ v[i];7 }8 c[row] = result ;9 barrier (CLK_LOCAL_MEM_FENCE);

10 }

Listing 4.11: Hilfsfunktion zur Matrix-Vektor-Multiplikation

Der Funktion werden die Adressen der Speicherbereiche übergeben, auf denengearbeitet wird. Hier wird wieder die Möglichkeit der parallelen Berechnungvon Matrixoperationen ausgenutzt. Wie in Listing 4.11 gezeigt, wird zunächsteine Variable result im Private Memory angelegt und mit null initialisiert.In einer Schleife wird die Summe aus den Produkten einer Spalte in resultzwischengespeichert und anschließend in der entsprechenden Spalte im Ergeb-nisvektor ~c abgelegt. Hier wird ebenfalls mittels des barrier()-Aufrufs sicher-gestellt, dass alle Ergebnisse berechnet sind, bevor mit der Ausführung desCodes fortgefahren wird, um Dateninkonsistenz zu vermeiden.Neben der in Listing 4.11 gezeigten Funktion existieren vier weitere Hilfsfunk-tionen für Matrix-Vektor-, beziehungsweise Matrix-Matrix-Operationen, die inTabelle 4.1 angegeben sind. Alle arbeiten nach dem zuvor beschriebenen Prin-zip. Der Unterschied zu der zuvor erläuterten Funktion besteht darin, dasseinige davon eine der Eingabevariablen überschreiben um Rechenoperationenzu sparen, falls diese nicht weiter benötigt wird.

c© Maik Dankel 44


4. ImplementierungNach Berechnung des Wertes pk|k−1 erfolgt ein erneuter Aufruf einer Unter-funktion, die die Kovarianzmatrix Ck|k−1 mittels Gleichung 3.2 berechnet. Dain diesem Fall Q = 0 gelten soll, ergibt sich daraus die vereinfachte Form

Ck|k−1 = FkCk−1|k−1F Tk (4.1)

Die Berechnung von Gleichung 4.1 erfolgt in zwei Schritten. Den ersten Schritt,in dem zunächst Fk ∗ Ck−1|k−1 berechnet wird, zeigt Listing 4.12.

1 ...2 __local scalar_t temp[ORDER∗ORDER];3 //F_k ∗ C_kk_14 for ( int i = 0; i < ORDER; i++) {5 result += jacobian[ORDER ∗ row + i] ∗ C_k_1k_1[column + ORDER ∗ i];6 }7 temp[ORDER ∗ row + column] = result;8 barrier (CLK_LOCAL_MEM_FENCE);9 ...

Listing 4.12: Erster Teil der Berechnung von Ck|k−1

Es wird ein temporäres Array temp angelegt, in dem das Ergebnis dieser Rech-nung zwischengespeichert wird. Die Schleife in Zeile 4 bis 6 führt die Ma-trixmultiplikation durch. Es wird für jede Stelle der Matrix die Summe derProdukte in der jeweiligen Zeile und Spalte gebildet und anschließend in derVariable temp gespeichert.Im Anschluss daran wird das Ergebnis dieses ersten Schrittes mit der transpo-nierten Jacobi-Matrix multipliziert, um das Endergebnis zu bilden.

1 ...2 //(F_k ∗ C_kk_1) ∗ F_k^T3 result = 0;4 for ( int i = 0; i < ORDER; i++) {5 result += temp[ORDER ∗ row + i] ∗ jacobian[ORDER ∗ column + i];6 }7 C_kk_1[ORDER ∗ row + column] = result;8 }

Listing 4.13: Zweiter Teil der Berechnung von Ck|k−1

Der zweite Teil dieser Funktion multipliziert das Zwischenergebnis des erstenTeils mit der transponierten Jacobi-Matrix. Mittels Manipulation der Indizes

c© Maik Dankel 45


4. Implementierungwird erreicht, dass keine weiteren Rechenoperationen nötig sind, um die Trans-ponierte zu erhalten. Diese einfache Anpassung erlaubt ohne weitere Änderungder Daten, und damit ohne höheren Reche

Implementierung eines GPU-beschleunigten Kalman-Filters ... · Implementierung eines...

Documents

Transcript of Implementierung eines GPU-beschleunigten Kalman-Filters ... · Implementierung eines...