GPU ProgrammierungEntwicklung paralleler Algorithmen mittels GPU

Hardwarebeschleunigung

Masterprojekt

Autoren: Philipp Schoppe, Matthias TöppeFachbereich: Elektrotechnik und Informatik

Betreuer: Prof. Dr. rer. nat. Nikolaus WulffDatum: 16. Februar 2014

AbstractIm Rahmen des Masterprojekts im Studiengang Master Informationstechnik undzweier Masterarbeiten ist ein Forschungsteam zum Thema GPU Programmierungzu begründen. Das initiale Team rund um Professor Dr. rer. nat. Nikolaus Wulffbesteht aus den Projektteilnehmern Philipp Schoppe und Matthias Töppe sowieMaik Dankel und Rene Böing, die durch ihre Forschungserkenntnisse, die sie wäh-rend der Erstellung ihrer Masterarbeiten und Promotion erlangen, zum Projektbeitragen. In der zweisemestrigen Projektphase ist die Infrastruktur des Laborsfür Informatik für GPU Programmierung aufzubauen und ein Einstieg in die Ent-wicklung von parallelen Algorithmen für Grafikkarten zu vollziehen.

i

Inhaltsverzeichnis

Abbildungsverzeichnis iv

1 Einleitung 11.1 GPU Computing . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 ATLAS Spurrekonstruktion per Kalman-Filter . . . . . . . . . . . . 2

1.2.1 Der ATLAS Detektor . . . . . . . . . . . . . . . . . . . . . . 21.2.2 Kalman-Filter . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.3 Infrastruktur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3.1 Projektserver . . . . . . . . . . . . . . . . . . . . . . . . . . 61.3.2 Laborrechner . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2 Kalman Filter Implementation 72.1 Datenextraktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 CUDA 123.1 Hardwaremodell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123.2 CUDA Programming Model . . . . . . . . . . . . . . . . . . . . . . 133.3 Speicherhierarchie . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153.4 Entwicklung paralleler Algorithmen . . . . . . . . . . . . . . . . . . 15

3.4.1 Vorgehensweise . . . . . . . . . . . . . . . . . . . . . . . . . 153.4.2 CUDA Best Practices . . . . . . . . . . . . . . . . . . . . . . 16

4 SkelCL 174.1 Skeletons . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174.2 Vector Datentyp . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

4.2.1 Speichermanagement . . . . . . . . . . . . . . . . . . . . . . 204.2.2 Distribution auf verschiedene Devices . . . . . . . . . . . . . 21

5 Lösen linearer Gleichungssysteme mit SkelCL 235.1 Gaußsches Eliminationsverfahren . . . . . . . . . . . . . . . . . . . 23

5.1.1 LU Dekomposition . . . . . . . . . . . . . . . . . . . . . . . 255.1.1.1 LU-Zerlegung . . . . . . . . . . . . . . . . . . . . . 255.1.1.2 Vorwärtseinsetzen . . . . . . . . . . . . . . . . . . 26

ii

Inhaltsverzeichnis

5.1.1.3 Rückwärtseinsetzen . . . . . . . . . . . . . . . . . . 265.2 Umsetzung mit SkelCL . . . . . . . . . . . . . . . . . . . . . . . . . 27

5.2.1 Vorwärtselimination/LU-Zerlegung . . . . . . . . . . . . . . 275.2.2 Vorwärtseinsetzen/Rückwärtseinsetzen . . . . . . . . . . . . 305.2.3 Nachiteration . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5.3 Ergebnisse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

6 Fazit und Ausblick 35

Literaturverzeichnis 36

iii

Abbildungsverzeichnis

1.1 Die GPU widmet der Datenverarbeitung mehr Transistoren. [1] . . 21.2 ATLAS Detektor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Messpunkte, Teilchenspur und Kalman-Filter . . . . . . . . . . . . . 41.4 Ablauf des Kalman Filter Verfahrens. . . . . . . . . . . . . . . . . . 4

3.1 CUDA Grid. [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133.2 Speicherhierarchie. . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

4.1 SkelCL - Architektur . . . . . . . . . . . . . . . . . . . . . . . . . . 174.2 Einfaches Beispiel eines SkelCL Skeletons [7] . . . . . . . . . . . . 204.3 Distributionsmethoden in SkelCL [7] . . . . . . . . . . . . . . . . . 22

5.1 Unabhängige Daten bei der Vorwärtselimination . . . . . . . . . . . 285.2 Unabhängige Daten bei dem Rückwärtseinsetzen . . . . . . . . . . . 30

iv

1 EinleitungDer Leistungszuwachs moderner Rechensysteme kann nicht mehr primär durchhöhere Taktraten erzielt werden. Infolgedessen rückt die Parallelisierung von Re-cheneinheiten in den Fokus der Architekturentwicklungen. Um dies auszunutzen,besteht die Notwendigkeit vorhandene Algorithmen zu parallelisieren oder neueBerechnungsverfahren zu entwerfen, die verteilte Arbeitsweisen ermöglichen.

Gängige parallele Recheneinheiten sind Multiprozessorsysteme (mit gemeinsamenoder verteilten Speicher), FPGAs (Field Programmable Gate Array) und GPU(Graphics Processing Unit) Systeme. Im Rahmen des Masterprojekts werden dieEigenschaften von GPU Systemen betrachtet. Abschnitt 1.1 gibt zunächst eineÜbersicht über das sogenannte GPU Computing, bevor Kapitel 3 und Kapitel 4 de-taillierter auf zwei Frameworks zur GPU Programmierung eingehen. Abschnitt 1.2stellt ein Forschungsthema im Bereich der Hochenergiephysik vor, bei dem derEinsatz von GPUs einen Geschwindigkeitsvorteil bringen soll. Dieses Thema warder Hauptbestandteil der beiden Masterarbeiten und wurde von dem Masterpro-jekt unterstützend begleitet. Danach folgt ein Abschnitt über die eingerichteteInfrastruktur für das Projekt.

1.1 GPU ComputingGPU Computing bezeichnet nach NVIDIA das durch Grafikprozessoren beschleu-nigte Rechnen. Der Grafikprozessor (GPU) kooperiert dabei mit der CPU undkann rechenintensive Anwendungen beschleunigen [2]. Grafikprozessoren wurdenfür Berechnungen in der Computergrafik entwickelt, in der häufig der gleiche Algo-rithmus auf sehr viele Pixelbereiche angewendet wird. Durch diesen Einsatzbereichbedingt, ist ihr Aufbau hoch parallel strukturiert: Es werden gegenüber einer CPUmehr Transistoren für die Datenverarbeitung eingebaut, wohingegen CPUs denDesignfokus auf einen guten Programmablauf und Caching legen (Abbildung 1.1).Moderne GPUs können über 2000 Recheneinheiten beinhalten. Sie sind optimiertfür Fließkommaoperationen und arbeiten seperat und unabhängig von der CPU.In einigen Forschungs- und Industriezweigen werden GPUs nicht mehr nur fürAufgaben in der Computergrafik sondern für allgemeine Rechenzwecke genutzt.Klassische Gebiete hierfür sind aufwändige Simulationen und Matrizenoperation.

1

1 Einleitung

Abbildung 1.1: Die GPU widmet der Datenverarbeitung mehr Transistoren. [1]

Grafikkartenhersteller ordnen diesem Trend das Schlagwort GPGPU (General Pur-pose Computation on Graphics Processing Unit) zu.

Es existieren verschiedene Frameworks zur Programmierung von GPUs. Die beidenbekanntesten sind das Open-Source PaketOpenCL und das proprietäre CUDA Fra-mework. OpenCL ist eine Erweiterung von C99 und drückt parallele Programmab-läufe aus. Eine Besonderheit an OpenCL ist, dass es für Produkte verschiedenerHersteller genutzt werden kann und auch auf Multiprozessorsystemen funktioniert.CUDA ist eine C-Erweiterung, die jedoch nur für Nvidia Grafikkarten vorgesehenist. Andere Ansätze verfolgen das Ziel ein Programm mit Pragma-Direktiven ver-sehen zu können, anhand deren der Compiler automatisch parallelen Code erzeugt.Beispiele hierfür sind OpenAAC und OpenMP. Eine auf OpenCL basierte SkeletonBibliothek ist SkelCL (siehe Kapitel 4). Sie soll die Grafikkartenprogrammierungdurch vordefinierte Templates vereinfachen.

1.2 ATLAS Spurrekonstruktion per Kalman-Filter1.2.1 Der ATLAS DetektorAm CERN (Europäische Organisation für Kernforschung) ist der LHC (Large Ha-dron Collider) beheimatet. Er ist der größte Teilchenbeschleuniger der Welt unddient zur Erforschung des Aufbaus von Materie. Die Teilchen innerhalb des LHCswerden auf nahezu Lichtgeschwindigkeit beschleunigt und zur Kollision gebracht.Dieser Moment wird von verschiedenen Detektoren festgehalten und analysiert. Ei-ner der Detektoren trägt den Namen ATLAS. Er untersucht Proton-Proton Kolli-sionen. Bei einer Kollision entstehen neue Teilchen die durch den Impuls der Kolli-sion verschiedene Detektorschichten passieren. Diese Spuren sind bei den Analysenvon wesentlicher Bedeutung. Grafik 1.2 skizziert ein solches Szenario. Da nur ein-zelne Messpunkte bekannt sind und keine komplette Spur, werden verschiedeneAlgorithmen angewandt, die aus den Sensordaten die entsprechenden Spuren ge-

2

1 Einleitung

Abbildung 1.2: An computer generated image representing how ATLAS detectsparticles. [6]

nerieren sollen. Ein Algorithmus, der zur Messfehlerminimierung eingesetzt wird,ist der im Folgenden vorgestellte Kalman-Filter.

1.2.2 Kalman-FilterDer Kalman-Filter ist ein von Rudolf E. Kalman vorgestellter Filter für dynami-sche Systeme. Als dynamisches System bezeichnet man das mathematische Modelleines zeitabhängigen Prozesses. Der Filter dient zur Analyse dieser dynamischenSysteme mithilfe einer Zustandsraummodellierung: Ein System wird insbesonde-re durch Zustände und zugehörige Zustandsübergänge modelliert [5]. Die Analyseermöglicht aus fehlerbehafteten Beobachtungen, welche aufgrund von Rauschenoder Ungenauigkeiten von Messgeräten auftreten, Rückschlüsse auf den tatsäch-lichen physikalischen Prozess zu ziehen. Anwendung findet dieses Verfahren zumBeispiel in der Kommunikationstechnik und Signalverarbeitung, in der Navigationoder bei der Wettervorhersage.

Ein weiteres Anwendungsgebiet liegt im Track- und Vertexfitting: In seinem Auf-satz „Application of Kalman filtering to track and vertex fitting“ [4] beschreibtFrühwirth den Zusammenhang zwischen Kalman-Filter und Trackfitting: Das Track-fitting dient dem Abschätzen der Spurparameter; der Kalman-Filter der Analyse(linearer) dynamischer Systeme. Betrachtet man eine Teilchenspur im (Detektor-)

3

1 Einleitung

Raum als dynamisches System, so kann das Filterverfahren als Fittingverfahrengenutzt werden, das die Messgenauigkeit nachträglich verbessert. Dies ist leichtmöglich, indem der Zustandsvektor des Systems als 5-dimensionaler Vektor p, derjeden Punkt der Spur eindeutig beschreibt, aufgefasst wird. Abbildung 1.4 skiz-ziert diese Idee:

Abbildung 1.3: Messpunkte der Detektorlagen, Teilchenspur und Filterung durchden Kalman-Filter. [3]

p könnte hierbei etwa durch eine Funktion mit einer passenden Koordinate z gebil-det werden, z.B.: p = p(z). In der Praxis werden die Schnittpunkte der Teilchen-spur mit den Detektorlagen als Koordinaten verwendet. Das Verfahren genügt einerrekursiven Prädiktor-Korrektor Struktur: Zunächst wird der zukünftige Systemzu-stand extrapoliert; in einem weiteren Schritt wird unter Betrachtung aller bishe-rigen Messungen diese Schätzung korrigiert und der aktuelle Zustand geschätzt.Für alle weiteren Messpunkte wird das Verfahren wiederholt. Das Verfahren isthinsichtlich der mittleren quadratischen Abweichung des Schätzers optimal [4].

Abbildung 1.4: Ablauf des Kalman Filter Verfahrens. Time update führt eine apriori Schätzung des Systemzustand durch. Das measurement up-date korrigiert dies durch eine a posteriori Schätzung. [9]

Die genauen mathematischen Formeln sind die folgenden:

4

1 Einleitung

Prädiktion

pk|k−1 = Fk−1pk−1|k−1 (1.1)

Extrapolation des Zustandvektors für die k-te Messung, mit dem vorausgehendenZustandsvektor pk−1|k−1 und der Jakobimatrix F.

Ck|k−1 = FkCk−1|k−1FTk + Qk (1.2)

Extrapolation der Kovarianzmatrix, mit der Jakobimatrix F, der vorausgehendenKovarianzmatrix Ck−1|k−1 und einem prozessbedingten Rauschterm Qk.

Korrektur (Gain Formalismus)

Die Extrapolationen aus der Prädiktionsphase werden in diesem Schritt unter Ein-fluss des aktuellen Messwerts korrigiert.

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

k )−1 (1.3)

Kk die Kalman gain Matrix, gibt an wie stark die Korrektur den neuen Zustandbeeinflusst. Mit der Transformationsmatrix Hk und dem Rauschterm Vk.

pk|k = pk|k−1 + Kk(mk −Hkkk|k−1) (1.4)

A posteriori Schätzung des Zustandsvektors unter Berücksichtigung der Kalmangain Matrix.

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

Korrektur der Kovarianzmatrix. Mit der Identitätsmatrix I.

Als weiterer Schritt des Filters kann das sogenannte Smoothing (Glättung) ange-wendet werden. Wurde eine Messreihe gefiltert kann abschließend für jeden Mess-punkt eine Zustandsschätzung durchgeführt werden, bei der die Informationen dergesamten Messreihe mit einfließen. Dieser Schritt verbessert die Schätzung derZustandsvektoren offensichtlich: Bei sehr „frühen“ Schätzungen einer Messreihefließen nur wenige Informationen vorausgehender Messpunkte in die Schätzungein - beim Smoothing hingegen alle Messpunkte. Die Implementationen der Mas-terarbeiten und des Masterprojekts realisieren dies, indem die Messreihe zunächstwie gewohnt vorwärts (vom ersten bis zum letzten Messpunkt) und in einem zwei-ten Schritt rückwärts (vom letzten bis zum ersten Messpunkt) gefiltert wird. DieZustandsschätzungen aus Hin- und Rückweg werden gewichtet zusammengeführt,sodass eine neue Zustandsschätzung entsteht. In dieser Schätzung fließen die Infor-mationen aller gegebenen Messpunkte ein. Im folgenden werden die Ergebnisse des

5

1 Einleitung

Hinwegs mit einem hochgestellten f (forward) bzw. einem hochgestellten b (back-ward) für den Rückweg gekennzeichnet. Daraus ergeben sich folgende Formeln:

Smoothing

K′k = Cfk|k(Cf

k|k + Cbk|k)−1 (1.6)

Formel zur Gewichtung der einzelnen Messwerte, bei die Kovarianzmatrizen ausHin- und Rückweg benutzt werden.

p′k = pfk|k + K′k(pb

k|k−1 − pfk|k) (1.7)

Zur Ermittlung des geglätteten Zustands wird das im Hinweg korrigierte pk|k mitdem im Rückweg extrapolierten pk|k−1 verbunden. Dies wirkt einer doppelten Ge-wichtung des korrigierten Werts entgegen.

1.3 Infrastruktur1.3.1 ProjektserverFür das Projekt wurde ein Server eingerichtet, der zum Sammeln und Festhal-ten von projektbezogenen Daten genutzt wird. Mithilfe der Software MediaWikiwurde ein Wiki installiert, das unter https://www.atlas.fh-muenster.deaufgerufen werden kann. Es ist jedoch nur für Projektteilnehmer mit eingerich-teten Benutzerkonten einsehbar. Ein SVN Repository ist unter https://www.atlas.fh-muenster.de/cuda erreichbar. Das cuda Projekt beinhaltet jeweilseine CUDA, OpenCL und CPU basierte Version.

1.3.2 LaborrechnerAuf den Arbeitsplätzen im Labor für Informatik wurden folgende Softwarepaketeinstalliert, die im Rahmen des Projekts benötigt wurden:• NVIDIA Grafikkartentreiber und CUDA Framework

• OpenCL

• SkelCL

• GNU Scientific Library

• CMake

• Subversion

• ROOT - A Data Analysis Framework

6

2 Kalman Filter ImplementationIm Rahmen des Masterprojekts wurde keine GPU Implementation des Kalman-Filters entwickelt. Als Unterstützung für die beiden Masterarbeiten wurde jedocheine CMake Buildumgebung, die Verarbeitungskette der zugrundeliegenden Datenund eine serielle Kalman-Filter Implementation zum Ergebnisvergleich program-miert. Der Zugriff auf die Daten erfolgt mit Hilfe des ROOT Frameworks. DieBerechnungen für den Kalman-Filter werden mit der GNU GSL Bibliothek voll-zogen.

2.1 DatenextraktionDie Daten der Events im ATLAS Detektor werden nach einigen Vorverarbeitungs-stufen in einem TTree Objekt des ROOT Frameworks gespeichert und serialisiert,sodass sie in eine ROOT Datei gespeichert werden können. Das ROOT Frame-work (http://root.cern.ch/) ist ein Framework zur Datenanalyse, das eineFülle von Algorithmen und Visualisierungsmöglichkeiten bereitstellt. Es ist sehrgut geeignet extrem große Datenmengen zu verwalten und zu kombinieren. DieEventdaten beinhalten durch Mustererkennung gefundene und vorgefittete Spu-ren auf denen der Kalman-Filter angewendet wird. Ein Event oder Ereignis be-schreibt eine Teilchenkollision. Die Zerfallsspuren innerhalb eines Ereignis werdenin Form von Koordinaten der Detektorlagen zur Verfügung gestellt. Ein TTreeObjekt besteht aus einer Liste von TBranch Objekten, die vom ROOT Frame-work adressiert werden können. Sie erlauben den Zugriff auf die Daten. Die einzel-nen Messpunkte der Detektorlagen sind in Form eines std::vector<float>Objekts abgespeichert. Sie bilden eine Spur. Spuren werden wiederum in einenVektor gekapselt; zusammen bilden sie ein Ereignis. Eine Branch besteht nun ausmehreren Entries, die jeweils ein Ereignis repräsentieren. Innerhalb eines Entriessteht folglich ein Objekt vom Typ std::vector<std::vector<float>> -Jedes Ereignis besteht aus diversen Spuren mit diversen Messpunkten. Eine wei-tere Komplexität besteht darin, dass die Messpunkte nicht als mehrdimensionaleDatentypen vorliegen, sondern für jede Komponente ein Branch existiert. So gibtes die Branches pseudoTrk_allhit_locX und pseudoTrk_allhit_locY. Damit nunall diese Daten in leichter verwaltbaren Strukturen verwendet werden können, wur-de die EventReader Klasse implementiert. Sie realisiert den Datenzugriff auf die

7

2 Kalman Filter Implementation

ROOT Datei und extrahiert die Daten in folgende Strukturen:

Listing 2.1: Strukturübersicht.struct TrackInfoStruct {

scalar_t d0;scalar_t z0;scalar_t phi;scalar_t theta;scalar_t qoverp;

};

struct TrackHitStruct {scalar_t normal[ORDER];scalar_t ref[ORDER];scalar_t err_locX;scalar_t err_locY;scalar_t cov_locXY;scalar_t jacobi[ORDER * ORDER];scalar_t jacobiInverse[ORDER * ORDER];char is2Dim;int detType;int bec;

};

typedef struct TrackHitStruct TrackHit_t;typedef struct TrackInfoStruct TrackInfo_t;

typedef std::vector<TrackHit_t> TrackData_t;

struct TrackStruct {TrackData_t track;TrackInfo_t info;TrackInfo_t truthTrackInfo;

};

typedef TrackStruct Track_t;

typedef std::vector<Track_t> KF_Event_t;

Ein Ereignis (KF_Event_t) besteht aus einer Spur track (TrackData_t) und In-formationen info (TrackInfo_t) über den Startpunkt der Spur. Falls die Datenaus einer Simulation stammen werden zusätzlich Informationen über den echtenStartpunkt aus der Simulation in das Feld truthTrackInfo hinterlegt.

2.2 ImplementationDie serielle Implementierung besteht aus einer Klasse KalmanFilterSerial,die jeweils eine Spur vom Typ Track_t einliest und den Filter auf ihr anwendet.Eine detaillierte Ausführung der gesamten Implementation wäre an dieser Stellezu umfangreich. Stattdessen werden lediglich ein paar Besonderheiten aufgezeigt.Beispielhaft sei hier die Methode zur Berechnung des Kalman-Gains vorgestellt:

8

2 Kalman Filter Implementation

Listing 2.2: Berechnung des Kalman-Gains.// Compute Kalman gain using the formula:// K_(k) = C_(k|k-1) H_(k)^T (V_(k) + H_(k) C_(k|k-1) H_(k)^T)^-1int KalmanFilterSerial::correctGain(GSL_MATRIX *dest, GSL_MATRIX *C_k1,

GSL_MATRIX *H_k, GSL_MATRIX *V_k) const {

if (!dest || !C_k1 || !H_k || !V_k) {PRINT_WARN("Parameter is NULL!");return -1;

}

#if DBG_LVL > 2std::cout << "correctGain:" << std::endl;std::cout << "H_k: ";printMatrix (H_k);std::cout << "C_k: ";printMatrix (C_k1);std::cout << "V_k: ";printMatrix (V_k);

#endif

GSL_MATRIX *tmp, *tmp2, *tmp3, *inverse;

tmp = GSL_MATRIX_CALLOC (H_k->size1, H_k->size2);tmp2 = GSL_MATRIX_CALLOC (H_k->size1, H_k->size1);tmp3 = GSL_MATRIX_CALLOC (H_k->size2, H_k->size1);inverse = GSL_MATRIX_CALLOC (H_k->size1, H_k->size1);

GSL_BLAS_XGEMM(CblasNoTrans, CblasNoTrans, 1.0, H_k, C_k1, 0.0, tmp);GSL_BLAS_XGEMM(CblasNoTrans, CblasTrans, 1.0, tmp, H_k, 0.0, tmp2);

GSL_MATRIX_ADD(tmp2, V_k);

GSL_MATRIX_MEMCPY(inverse, tmp2);

matrixInverseSymmetric(inverse);

GSL_BLAS_XGEMM(CblasTrans, CblasNoTrans, 1.0, H_k, inverse, 0.0, tmp3);GSL_BLAS_XGEMM(CblasNoTrans, CblasNoTrans, 1.0, C_k1, tmp3, 0.0, dest);

GSL_MATRIX_FREE(inverse);GSL_MATRIX_FREE(tmp);GSL_MATRIX_FREE(tmp2);GSL_MATRIX_FREE(tmp3);

#if DBG_LVL > 2std::cout << "maxima correct Gain:" << std::endl;std::cout << "C_k . transpose (H_k) . invert ((V_k + H_k . C_k . transpose (H_k)));";std::cout << std::endl;std::cout << "res: ";printMatrix (dest);

#endif

return 0;}

An dieser Stelle ist zu sehen, dass die BLAS Routinen der GNU GSL Biblio-thek durch Makros gekapselt werden. Sie ermöglichen einen einfachen Wechselzwischen float und double Genauigkeiten. Normalerweise sieht die Biblio-thek unterschiedliche Prototypen der Routinen vor. Statt gsl_blas_dgemm undgsl_blas_sgemm gibt es nur ein Makro GSL_BLAS_XGEMM. Ein weiteres Merk-mal besteht in den Debugausgaben: Sie sind so strukturiert, dass alle Berechnun-gen durch simples kopieren mit dem Computeralgebrasystem maxima überprüftwerden können. Der Aufruf der Methode matrixInverseSymmetric deutet

9

2 Kalman Filter Implementation

bereits an, dass zu verarbeitende Matrix immer symmetrisch ist. Daher kann dieCholesky-Zerlegung zur Invertierung genutzt werden. Anfangs wurde eine normaleInvertierung durchgeführt, die Cholesky-Zerlegung ist jedoch bedeutsam schneller.Listing 2.3 und 2.4 zeigen beide Varianten. Die allgemeine Matrixinvertierung:

Listing 2.3: Allgemeine Matrixinvertierung.// Calculate the inverse of a matrix src.inline void matrixInverse (GSL_MATRIX *src) {

gsl_matrix * LU = gsl_matrix_calloc (src->size1, src->size2);gsl_matrix * tmp = gsl_matrix_calloc (src->size1, src->size2);

for (size_t i = 0; i < src->size1; ++i) {for (size_t j = 0; j < src->size2; ++j)

gsl_matrix_set (LU, i, j, (double) GSL_MATRIX_GET (src, i, j));}

int s;gsl_permutation * perm = gsl_permutation_calloc (LU->size1);gsl_linalg_LU_decomp (LU, perm, &s);

// LU decomp singular?for (size_t i = 0; i < LU->size1; i++) {

if (gsl_matrix_get (LU, i, i) == 0) {

#if DBG_LVL > 2std::cout << "Matrix is singular! Using pseudoinverse matrix." << std::endl;

#endifgsl_matrix_free (LU);gsl_permutation_free (perm);gsl_matrix_free (tmp);

pseudoInverse (src);return;

}}

gsl_linalg_LU_invert (LU, perm, tmp);

for (size_t i = 0; i < src->size1; ++i) {for (size_t j = 0; j < src->size2; ++j)

GSL_MATRIX_SET (src, i, j, (scalar_t) gsl_matrix_get (tmp, i, j));}

gsl_matrix_free (LU);gsl_matrix_free (tmp);gsl_permutation_free (perm);

return;}

Im Fall der allgemeinen Invertierung wird auch überprüft, ob die Matrix singu-lär ist. In diesem Fall wird eine Pseudoinverse berechnet. Ohne diese Überprü-fung würde die standardmäßige Fehlerbehandlung der GSL Bibliothek greifen unddas Programm beenden. Während bei der allgemeinen Invertierung noch mehre-re Schritte im Zusammenspiel mit einer LR-Zerlegung vollzogen werden müssen,wird die schnellere Invertierung symmetrischer Matrizen komplett von der GSLBibliothek unterstützt:

10

2 Kalman Filter Implementation

Listing 2.4: Matrixinvertierung per Cholesky-Zerlegung.inline void matrixInverseSymmetric(GSL_MATRIX *src) {

// Cholesky decomp only works with gsl_matrix (no float support)#if USE_DOUBLE

gsl_linalg_cholesky_decomp (src);gsl_linalg_cholesky_invert (src);

#elsegsl_matrix *tmp = gsl_matrix_calloc (src->size1, src->size2);

for (size_t i = 0; i < src->size1; ++i) {for (size_t j = 0; j < src->size2; ++j)

gsl_matrix_set (tmp, i, j, (double) GSL_MATRIX_GET (src, i, j));}

gsl_linalg_cholesky_decomp (tmp);gsl_linalg_cholesky_invert (tmp);

for (size_t i = 0; i < tmp->size1; ++i) {for (size_t j = 0; j < tmp->size2; ++j)

GSL_MATRIX_SET (src, i, j, (scalar_t) gsl_matrix_get (tmp, i, j));}

gsl_matrix_free (tmp);#endif

return;}

11

3 CUDACUDA steht für eine GPGPU Plattform für NVIDIA Grafikkarten. Sie umfasstsowohl Treiber für die Grafikkarten als auch eine Runtime API, die die Program-mierung für GPUs erlaubt. Indem sie die Programmiersprache C um nur weni-ge zusätzliche Konstrukte erweitert, wird der Einstieg für Entwickler erleichtert.Nachdem zunächst das Hardwaremodell (siehe Abschnitt 3.1) vorgestellt wird, wer-den in Abschnitt 3.2 diese Konstrukte und ihre Logik erläutert. Abschließend wirdvorgestellt wie sich Entwickler ohne vorhandenes Spezialwissen der GPGPU Weltund parallelen Algorithmen nähern können.

3.1 HardwaremodellEine NVIDIA GPU besteht aus mehreren verschiedenen Hardwarekomponenten.Neue Grafikkarten sind über ein PCI Express 3.0 Interface mit dem Hostsystemverbunden. Das Hostsystem kann somit Daten in den GPU Speicher kopieren undden Start von sogenannten Kernelfunktionen veranlassen. Dies sind von Entwick-lern programmierte Funktionen die von dem Scheduler an ein skalierbares Feld vonmultithreaded Streaming Multiprocessors (SMs) verteilt werden. Sie sind gewis-sermaßen mit CPUs vergleichbar: Sie besitzen ebenfalls Register, einen Instruk-tionsspeicher, diverse Rechenkerne und einen Shared Memory Bereich über dendie Threads Daten austauschen können. Ein SM ist so entworfen, dass er hunder-te von Threads parallel ausführen kann. Dies wird durch die sogenannte SIMT(Single-Instruction, Multiple Thread) Architektur erreicht, bei der ein Schedulerjeweils 32 Threads gruppiert. Die Threads werden zusammen erstellt, kontrolliertund können in einem Schritt die gleiche Anweisung ausführen. Diese Gruppierungwird Warp genannt. Die einzelnen Threads starten zusammen an der gleichenProgrammspeicheradresse, verfügen aber über eigene Programmzähler und Regis-ter. Im Optimalfall führt jeder Thread eines Warps die gleiche Instruktion aus,damit die höchstmögliche Performance erreicht werden kann. Sollten Programm-sprünge aufgrund von datenabhängigen Programmanweisungen auftreten, werdendie Anweisungen der Threads nicht mehr gleichzeitig ausgeführt, sondern unterPerformanceeinbußungen serialisiert. Dieses Verhalten ist von den Entwicklern zubeachten, um Effizient zu Programmieren. Die bereits erwähnten Kernel arbeitenauf logischen Datenblöcken, die der Scheduler dann auf die SMs verteilt. Die SMs

12

3 CUDA

teilen die Blöcke wiederum auf verschiedene Warps auf. Die Logik der Datenblöckewird im nächsten Abschnitt genauer erläutert.

3.2 CUDA Programming ModelDie logische Strukturierung von Daten entscheidet maßgeblich darüber, wie leichtAlgorithmen für abstrahierte Systeme erstellt werden können. Das CUDA Pro-gramming Model unterstützt die Abstraktion, indem es die SM Arrays als 3-dimensionales Feld repräsentiert. Die Threadgruppen werden in Blöcke eingeteilt,welche wiederum ein Raster, das Grid genannt wird, bilden: In einem Programm ist

Abbildung 3.1: CUDA Grid. [1]

es nicht zwingend die Logik des 3-dimensionalen Grids zu übernehmen. Je nach Pa-rametrisierung der Kernel kann auch auf einem (logischen) 1- oder 2-dimensionalenRaster gearbeitet werden. Im Folgendem wird die Syntax der Kernel erläutert undwie diese die einzelnen Threads verwalten. Die Gridstruktur sollte dadurch ver-ständlich werden. Im Gegensatz zu einer normalen C Funktion wird ein Kernel Nmal parallel von N verschiedenen Threads ausgeführt. Er wird durch das Schlüs-selwort __global__ ausgezeichnet. Die Anzahl der zu nutzenden Threads wirddurch die <<<...>>> Syntax spezifiziert. Listing 3.1 zeigt einen Kernel der ei-ne parallele Vektoraddition unter Ausnutzung von N Threads implementiert. Diesentspricht einer 1-dimensionalen Betrachtungsweise des Grids (es wird nur die ’x-Dimension’ verwendet). Jeder Thread der nun einen Kernel ausführt, besitzt eine

13

3 CUDA

eindeutige Identifikationsnummer. Aus Sicht eines Kernels ist sie in einer globalenVariablen threadIdx gespeichert. threadIdx ist vom Datentyp dim3, der einex, y und z Komponente beinhaltet. Je nachdem mit welchen Parametern ein Kernelgestartet wird, erlaubt dieser Datentyp eine leichte Zuordnung der Blockkompo-nenten zu Thread-IDs: Für einen 1-dimensionalen Block entspricht die Thread-IDder x Komponente; Für einen zweidimensionalen Block (Dx, Dy) ist die Thread-IDder Komponente (x,y) mittels (x + y Dx) zu bestimmen; Die ID der Komponente(x,y,z) eines Blocks (Dx, Dy, Dz) wird durch die Formel (x + y Dx + z Dx Dy) er-rechnet. Die Zuordnung von Thread-ID und der globalen Variablen ist ein CUDAinterner Mechanismus.

Listing 3.1: Vektoraddition mit N Threads. [1]// Kernel definition__global__ void VecAdd(float* A, float* B, float* C) {

int i = threadIdx.x;C[i] = A[i] + B[i];

}

int main() {...// Kernel invocation with N threadsVecAdd<<<1, N>>>(A, B, C);...

}

Für die Vektoraddition genügt die 1-dimensionale Betrachtungsweise. Für eineMatrix hingegen ist eine 2-dimensionale Betrachtungsweise sinnvoller:

Listing 3.2: Matrixaddition der NxN Matrizen A und B. [1]// Kernel definition__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]) {

int i = threadIdx.x;int j = threadIdx.y;C[i][j] = A[i][j] + B[i][j];

}

int main() {...// Kernel invocation with one block of N * N * 1 threadsint numBlocks = 1;dim3 threadsPerBlock(N, N);MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C);...

}

Bei der Entwicklung paralleler Algorithmen sollten die zugrunde liegenden Da-tenstrukturen so gewählt werden, dass sie mit dieser Programmierlogik leicht undübersichtlich adressierbar sind.

14

3 CUDA

3.3 SpeicherhierarchieIm allgemeinen besitzen NVIDIA Grafikkarten die folgende Speicherhierarchie:Jeder Thread hat einen eigenen Registerspeicher, auf den nur der einzelne Thread

Abbildung 3.2: Speicherhierarchie.

zugreifen kann. Alle Threads innerhalb eines Blocks haben die Möglichkeit übereinen Shared-Memory Bereich Daten untereinander auszutauschen. Zudem gibtes einen größeren globalen Speicherbereich, auf den jeder Thread in jedem BlockZugreifen kann. Der Speicherzugriff auf Register hat geringere Latenzzeiten als einZugriff auf den Shared-Memory. Zugriffe auf den Shared-Memory sind wiederumschneller als Zugriffe auf den globalen Speicher.

3.4 Entwicklung paralleler AlgorithmenDie Entwicklung paralleler Algorithmen ist ohne gute Vorkenntnisse der GPU Ar-chitektur nicht möglich. In Abschnitt 3.4.1 wird eine Vorgehensweise vorgestellt,wie seriell implementierte existierende Algorithmen mit Hilfe von CUDA beschleu-nigt werden können. Abschnitt 3.4.2 gibt einen Überblick über Optimierungsmög-lichkeiten für CUDA Programme.

3.4.1 VorgehensweiseEntwickler können mit Hilfe des Assess, Parallelize, Optimize, Deploy (APOD)Entwurfskreislaufs existierende Lösungen systematisch analysieren parallelisieren.Die Phasen werden der Reihe nach durchlaufen. Nach jedem Durchlauf kann derFokus auf immer feinere Details pro Zyklus gelegt werden. Das Ziel der einzelnenPhasen wird im Folgenden erläuter:

Assess

In der ersten Phase werden Programmabschnitte identifiziert, die den Großteil derGesamtlaufzeit ausmachen. Ist ein solcher Flaschenhals gefunden worden, muss

15

3 CUDA

der Entwickler die Grenzen der möglichen Optimierungen erfassen, indem er mitAmdahls und Gustafsons Gesetz untersucht wird. Sie geben Auskunft über dentheoretisch möglichen Geschwindigkeitsvorteil durch Parallelisierung.

Parallelize

Wurde festgestellt, dass ein identifizierter Flaschenhals durch Parallelisierung be-schleunigt werden kann, wird in dieser Phase der CUDA Code geschrieben. Viel-leicht stehen sogar schon fertige Methoden zur Verfügung, sodass lediglich einFunktionsaufruf geändert werden muss. Dies könnte etwas bei BLAS Routinen derFall sein. Je nach Programmdesign müssen eventuell aber auch zugrundeliegendeDatenmodelle angepasst werden. Es liegt am Entwickler den Programmieraufwandund möglichen Leistungsgewinn abzuwägen.

Optimize

Nach jeder Parallelize-Phase folgt eine Optimierungsphase. In dieser Phase wirdder CUDA Code auf potenzielle Optimierungsmöglichkeiten untersucht. Das Zu-sammenspiel zwischen bestehenden Programmteilen und neuen Kernelaufrufenmuss ebenfalls in Einklang gebracht werden. Zu den Optimierungen gehören sowohlAnpassungen am gesamten Programmfluss, wie Datentransfers oder Datenstruk-turen als auch feine Optimierungen wie die Optimierung von zusammenhängendenSequenzen von Fließkommaoperationen.

Deploy

Nachdem die ersten drei Phasen abgeschlossen sind, muss evaluiert werden, ob dieAnpassungen den zuvor bestimmten Erwartungen entsprechen. Der tatsächlicheund in der Assess-Phase bestimmte theoretische Geschwindigkeitsvorteil solltenverglichen werden. Kann die Anpassung als Erfolg gewertet werden, sollte sie nachbestanden Fach- und Qualitätstests so früh wie möglich produktiv eingesetzt wer-den.

3.4.2 CUDA Best PracticesIn der Optimierungsphase kann auf den CUDA Best Practices Guide (http://docs.nvidia.com/cuda/cuda-c-best-practices-guide/) zurückgegrif-fen werden. Er stellt eine Art Checkliste von Optimierungsmöglichkeiten dar, dievon jedem Entwickler beachtet werden sollte. Er gibt Hinweise zur Optimierungvon Datentransfers, optimale Speicherzugriffe auf den globalen und gemeinsamenSpeicherbereich sowie Hinweise zur Auslastung der Karte. Weniger ins Gewichtfallen Anweisungsoptimierungen, aber auch hierzu gibt es Tipps.

16

4 SkelCLDie an der Universität Münster entwickelte SkelCL Bibliothek ist eine High-LevelAPI basierend auf OpenCL und stellt vier vorimplementierte Berechnungs- undKommunikationsmuster, die sogenannten Skeletons, bereit. Durch diese Skeletonssoll die Entwicklung für GPU Systeme und speziell auch für Multi-GPU Systeme,also Systeme mit mehreren GPUs, erheblich vereinfacht und beschleunigt werden.Weiterhin stellt die SkelCL Bibliothek einen abstrakten Vector Datentypen sowieHigh-Level Mechanismen für die Distribution der Daten auf verschiedene GPUsbereit. Mit diesen Mechanismen und dem abstrakten Vector Datentyp werden dieDatentransfers zwischen dem Hauptspeicher und dem Speicher der einzelnen GPUsgekapselt und somit vor dem Entwickler verborgen. Architektonisch betrachtet bil-det die SkelCL Bibliothek somit eine eigene Schicht oberhalb von OpenCL, wie inAbbildung 4.1 dargestellt. [8] [7]

4.1 SkeletonsEin Skeleton ist im Grunde genommen ein Gerüst, dass eine BenutzerdefinierteFunktionen auf einer Menge von Daten anwenden kann. Ausgeführt werden dieFunktionen dann parallel auf einem oder mehrere OpenCL Devices. Dabei werdendie Details, wie z.B. mit den einzelnen Devices kommuniziert wird oder auch an-dere Technische Details, vor dem Benutzer verborgen.

Abbildung 4.1: SkelCL - Architektur

17

4 SkelCL

SkelCL bietet die vier Skeletons: map, zip, reduce und scan.

• Dasmap Skeleton wendet auf jedes Element eines Eingangsvektors [x1, ..., xn]die vom Benutzer definierte Funktion f an, so dass gilt:

map(f)([x1, ..., xn]) = [f(x1), ..., f(xn)] (4.1)

• Das zip Skeleton operiert auf zwei Eingangsvektoren ([x1, . . . , xn], [y1, . . . , yn])und wendet, einen vom Benutzer definierten, assoziativen binären Operator⊕ paarweise auf die Elemente der Eingangsvektoren an, so dass gilt:

zip(⊕)([x1, . . . , xn], [y1, . . . , yn]) = [(x1 ⊕ y1), ..., (xn ⊕ yn)] (4.2)

• Das reduce Skeleton errechnet einen skalaren Wert aus einem Eingangsvektor[x1, . . . , xn] indem es einen binären, vom Benutzer spezifizierten Operator �hintereinander auf alle Elemente des Vektors anwendet, so dass gilt:

reduce(�)([x1, ..., xn]) = x1 � x2 � ...� xn (4.3)

• Das scan Skeleton errechnet einen Vektor aus Zwischensummen indem eseinen assoziativen binären, vom Benutzer spezifizierten Operator � auf einenEingangsvektor [x1, . . . , xn] anwendet, so dass gilt:

scan(�)([x1, ..., xn]) = [x1, (x1�x2), (x1�x2�x3), ..., (x1�x2�...�xn)] (4.4)

Es ist ebenfalls möglich die einzelnen Skeletons mit weiteren Eingabeparameternzu versehen. Diese können entweder statische skalare Werte, die für jede Berech-nung identisch sind, oder aber weitere Eingabevektoren sein. In dem folgendenBeispiel in Abbildung 4.2 wird die Verwendung von zusätzlichen Parametern ver-deutlicht.

Alle vier Skeletons folgen dem Singe Instruction, Multiple Data (SIMP) Ausfüh-rungsmodell. Es wird also jeweils eine Funktion definiert, die auf multiplen Datenangewendet wird. Die Entwickler der Bibliothek haben sich für genau diese vierSkeletons entschieden, da diese in den Augen der Entwickler für eine breite Massean Applikationen hilfreich seien können.

Durch den Einsatz der SkelCL Skeletons müssen Entwickler nun mehr keine Low-Level Kernel in OpenCL implementieren, sondern sie können aus den Skeletonsein oder mehrere passende Skeletons auswählen und an die eigenen Bedürfnisse

18

4 SkelCL

anpassen. Die technischen Details der OpenCL Implementation wie z.B. das Allo-kieren von Speicherbereichen auf der GPU, deren Übergabe als Pointer, oder dasTransferieren der Daten von und zu der GPU werden von SkelCL realisiert undbleiben vor dem Entwickler verborgen. Durch dieses höhere Abstraktionsniveauwird nicht nur die Entwicklung erheblich beschleunigt, auch die Komplexität derAnwendung sinkt und somit sinkt auch die die Wahrscheinlichkeit, dass Fehler inder Anwendung implementiert werden.

Abbildung 4.2 zeigt wie einfach und schnell diese Skeletons angewendet und indi-viduell an die eigenen Bedürfnisse angepasst werden können. Dabei sind kaum biskeine technischen Details vom Entwickler zu beachten. Ein Entwickler kann so,mit sehr geringem Aufwand, die Rechenleistung seiner GPU bzw. seiner GPUs inseiner Anwendung nutzen.

In diesem Beispiel wird zunächst mit Zip<float> saxpy() ein neues zip Skeleton ineinfacher Genauigkeit (Single Precision) mit dem Namen saxpy erzeugt. Als Pa-rameter erhält dieser Aufruf die benutzerdefinierte und auszuführende Operationals String. Dieser String wird dann zur Laufzeit in Maschinencode für das entspre-chende Device übersetzt. In diesem Beispiel lautet die Funktionsdefinition:"float func(float x, float y, float a) return a*x+y; ". Die Funktion erwartet dreiParameter vom Typ float, hier sehen wir bereits, dass das zip Skeleton um einenweiteren Parameter a ergänzt wurde, denn das zip Skeleton operiert üblicherweiseauf zwei Eingangsvektoren. Der Parameter a ist in diesem Beispiel ein einfacherskalarer Wert. Die Operation ⊕ wird spezifiziert als: ax + y. Daraus folgt für dasErgebnis:zip(ax + y)([x1, . . . , xn], [y1, . . . , yn], a) = [(ax1 + y1), ..., (axn + yn)].Anschließend werden die Eingabevektoren X und Y befüllt und a ein skalarerWert zugewiesen. Abschließend wird durch die Zeile Y = saxpy(X, Y, a); die Aus-führung des Skeletons auf den vorhandenen GPUs angewiesen und die Ergebnisseder Berechnungen wiederum in dem Vektor Y abgelegt.

Wie zu erkennen ist, wurden keine Kenntnisse über die technischen Details derOpenCL Implementation benötigt und es ist mit sehr wenigen Zeilen Code mög-lich die vorhandenen GPUs für generische Berechnungen zu verwenden. Die SkelCLBibliothek generiert dabei zur Laufzeit, aus den vordefinierten Skeletons mit dengetroffenen Anpassungen, gültigen OpenCL Code und führt diesen transparent imHintergrund aus.[8] [7]

19

4 SkelCL

Abbildung 4.2: Einfaches Beispiel eines SkelCL Skeletons [7]

4.2 Vector DatentypDer von SkelCL implementierte Vector Datentyp stellt einen zusammenhängendenSpeicherbereich bereit, der sowohl von den einzelnen Devices als auch von dem Hostadressiert werden kann. In den Vectoren werden daher die zur Berechnung nötigenDaten gehalten und ebenfalls die Ergebnisse der Berechnungen abgelegt.

Eine weitere wichtige Funktion des Vectors ist die Verteilung (Distribution) derDaten und somit der Rechenaufgabe auf die verschiedenen Devices, sofern meh-rere vorhanden sind. So kann die vorhandene Performance optimal ausgenutztwerden.[8] [7]

4.2.1 SpeichermanagementIn OpenCL sind die Speicherbereiche des Hosts und der einzelnen Devices vonein-ander getrennt. Daher ist es notwendig die Daten explizit vom Host zu einem odermehreren Devices zu übertragen (Upload) oder umgekehrt (Download). Der Vectorkapselt dieses Speichermanagement und verbirgt die Details vor dem Entwicklerum die Entwicklung zu vereinfachen. Dazu verwaltet der Vector intern Pointerauf verschiedene Speicherbereiche der Devices und des Hosts. Sobald die Datendes Vectors entweder von dem Host oder von einem Device angefordert werden,veranlasst der Vector die entsprechende Übertragung der Daten und hält somit

20

4 SkelCL

die Daten zwischen den Devices und dem Host transparent im Hintergrund syn-chron. Daher ist der Vector in der Lage seine Daten sowohl dem Host als auchden einzelnen Devices zur Verfügung zu stellen. Der Datenaustausch erfolgt dabeilazy, sprich nur wenn es unbedingt notwendig ist. Wird beispielsweise der Ausgabe-Vector eines Skeletons als Eingabe-Vector eines zweiten Skeletons verwendet, dabeiallerdings nie vom Host gelesen oder verwendet, so verbleiben die Daten direkt aufdem Device und werden niemals zurück zum Host transferiert. Dadurch werdenkeine unnötigen Übertragungen ausgeführt die die Rechenleistung negativ beein-flussen könnten. [8] [7]

4.2.2 Distribution auf verschiedene DevicesBei Systemen mit mehreren verfügbaren Devices ist es sinnvoll, die Rechenaufgabeauf alle verfügbaren Devices aufzuteilen, da sich die Ausführungszeiten hierdurcherheblich verringern lassen. Daher realisiert der Vector Datentyp der SkelCL Bi-bliothek ebenfalls die Distribution der Daten und der Aufgaben auf die einzelnenDevices. Dazu implementiert der Vector die drei verschiedenen und frei wählbarenDistributionsmethoden: single, block und copy.

Abbildung 4.3 stellt die drei Distributionsmethoden grafisch dar. Wie dort zuerkennen ist, werden bei der Distributionsmethode single die gesamten Daten anein einzelnes Device gesendet. In der Regel ist dies das erste Device. Allerdingskann diese Vorbelegung auch überschrieben werden. Zu beachten ist allerdings,dass für alle Eingabe- und Ausgabevektoren eines Skeletons dasselbe Device alsZiel der Distribution gewählt werden muss. In dieser Methode wird also nur dieRechenleistung eines Devices verwendet.

Bei der Distributionsmethode block hingegen wird die Datenmenge in einzelneBlöcke geteilt und auf die vorhandenen Devices aufgeteilt. Jedes Device rechnetdann auf dem ihm zugeteilten Datenblock. Im Anschluss an die Berechnung wirddann aus den Teilergebnisblöcken das Gesamtergebnis zusammengesetzt. Bei die-ser Distributionsmethode wird eine optimale Laufzeit erzielt. Für alle Ein- undAusgabevektoren eines Skeletons, die block als Distributionsmethode verwenden,muss die gleiche Aufteilung in Blöcke gelten.

Bei der dritten und letzten Distributionsmethode copy erhält jedes Device eineneigenen vollständigen Datensatz, daraus folgt auch, dass alle Devices dieselbe Auf-gabe mit denselben Daten rechnen.

21

4 SkelCL

Abbildung 4.3: Distributionsmethoden in SkelCL [7]

Jedem Vector kann entweder explizit eine Distributionsmethode zugewiesen wer-den oder aber die Zuweisung einer Distributionsmethode erfolgt implizit durch dasSystem. Dazu setzt das System in Abhängigkeit von dem gewählten Skeleton dievermeintlich optimale Distributionsmethode. [8] [7]

22

5 Lösen linearer Gleichungssystememit SkelCL

Im zweiten Semester des Masterprojekts haben wir uns mit dem Lösen von linearenGleichungssystemen beschäftigt. Das Lösen von linearen Gleichungssystemen istein wichtiger Bestandteil sehr vieler Algorithmen, daher soll das Lösen der linearenGleichungssysteme soweit möglich parallelisiert werden und mittels der SkelCLBibliothek auf den leistungsfähigen GPUs implementiert werden.

5.1 Gaußsches EliminationsverfahrenDas gaußsche Eliminationsverfahren oder auch Gauß-Verfahren wurde von CarlFriedrich Gauß entwickelt und ist ein Algorithmus zum lösen linearer Gleichungs-systeme. Das Verfahren beruht darauf, dass elementare Umformungen zwar dasGleichungssystem ändern aber die Lösung dabei erhalten bleibt. Dies erlaubt esjedes eindeutig lösbare Gleichungssystem auf eine Stufenform zu bringen, bei derdurch sukzessive Elimination der Unbekannten, die Lösung mit geringem Aufwanderrechnet werden kann.

Ein lineares Gleichungssystem Ax = b mit n Gleichungen und m Unbekanntenx = (x1, ..., xm)T und b = (b1, ..., bn)T hat die Form:

a11x1 + a12x2 + ... + a1mxm = b1,a21x1 + a22x2 + ... + a2mxm = b2,

... + ... + . . . + ... = ...,an1x1 + an2x2 + ... + anmxm = bn.

(5.1)

oder:

a11 a12 ... a1m

a21 a22 ... a2m... ... . . . ...

an1 an2 ... anm

x1x2...

xm

=

b1b2...

bn

(5.2)

23

5 Lösen linearer Gleichungssysteme mit SkelCL

Für eine bessere Übersichtlichkeit können die Koeffizienten aij des Gleichungs-systems in der mit b erweiterten Koeffizientenmatrix dargestellt werden.

a11 a12 ... a1m b1a21 a22 ... a2m b2... ... . . . ... ...

an1 an2 ... anm bn

(5.3)

Gesucht wird eine Lösung für den Vektor x = (x1, ..., xm)T des lineares Gleichungs-system Ax = b.

Mit dem Gauß-Verfahren wird die gesuchte Lösung für x in zwei Schritten be-rechnet:

1. Vorwärtselimination.

2. Rückwärtseinsetzen (Rücksubstitution).

Im ersten Schritt wird das Gleichungssystem auf die Stufenform gebracht, so dassin jeder Zeile mindestens eine Variable weniger auftritt als in der vorherigen Zeile.Man spricht auch davon, dass Variablen eliminiert werden. Erreicht werden kanndiese Elimination durch zwei Arten von elementaren Umformungen:

1. Zu einer Zeile wird ein vielfaches einer anderen Zeile addiert.

2. Zwei Zeilen werden miteinander vertauscht.

Diese elementaren Umformungen transformieren zwar das Gleichungssystem, er-halten allerdings die Lösungsmenge. Nach der Transformation soll die um b erwei-terte Koeffizientenmatrix folgende Form haben:

a′11 a′12 ... a′1m b′10 a′22 ... a′2m b′2... . . . . . . ... ...0 ... 0 a′nm b′n

(5.4)

Um diese Form zu erhalten iteriert der Algorithmus über die Spalten i = 1, . . . , m

24

5 Lösen linearer Gleichungssysteme mit SkelCL

und addiert zu jeder Zeile j = i + 1, . . . , n also der Zeilen unterhalb dem Diago-nalelement ein vielfaches der i-ten Zeile hinzu, so dass a′ij = 0 wird.

Im zweiten Schritt des Gauß-Verfahren, dem Rückwärtseinsetzen oder auch Rück-substitution genannt, wird schließlich die gesuchte Lösung für den Vektor x =(x1, ..., xm)T des lineares Gleichungssystems berechnet. Dazu wird beginnend aberder letzten Zeile die Unbekannten aufgelöst durch:

xi = 1a′ii

b′i −n∑

k=i+1a′ik · xk

(5.5)

5.1.1 LU DekompositionDie LU Dekomposition ist eine Variante des Gauß-Verfahrens welches häufig ver-wendet wird wenn quadratische Gleichungssysteme durch ein Computerprogrammgelöst werden sollen. Der Algorithmus arbeitet in drei Schritten:

1. LU-Zerlegung.

2. Vorwärtseinsetzen.

3. Rückwärtseinsetzen.

5.1.1.1 LU-Zerlegung

Zur Lösung des Gleichungssystems wird zunächst die Matrix A in das Produkteiner unteren Dreiecksmatrix L und eine oberen Dreiecksmatrix U zerlegt, so dassgilt:

A = LU (5.6)mit:

a11 a12 ... a1n

a21 a22 ... a2n... ... . . . ...

an1 an2 ... ann

=

l11 0 ... 0l21 l22

. . . ...... ... . . . 0

ln1 ln2 ... lnn

·

u11 u12 ... u1n

0 u22 ... u2n... . . . . . . ...0 . . . 0 unn

(5.7)

U entspricht dabei der oben erwähnten Stufenform des Gauß-Verfahrens. Dabei ist

25

5 Lösen linearer Gleichungssysteme mit SkelCL

anzumerken, dass hier nicht die um b erweiterte Koeffizientenmatrix auf die Stu-fenform gebracht wurde, sondern lediglich die Matrix A. In L werden die jeweiligenvielfachen der zu addierenden Zeilen, also die Umformungsschritte, hinterlegt. UmEindeutigkeit zu erreichen werden die Elemente Hauptdiagonalen von L gleich 1gesetzt.

5.1.1.2 Vorwärtseinsetzen

Da bei der LU Dekomposition der Lösungsvektor b = (b1, ..., bn)T nicht zusammenmit der Matrix A transformiert wurde, ist in diesem Verfahren ein zusätzlicherSchritt notwendig. Die benötigten Umformungsschritte um den Vektor b zu trans-formieren wurden in L gespeichert.Aus A = LU folgt:

Ly = b (5.8)und

Ux = y (5.9)Wir müssen also zunächst das Gleichungssystem Ly = b lösen. Da auch L eineDreiecksmatrix ist, kann die Lösung ermittelt werden durch:

yi = 1lii

(bi −

i−1∑k=1

lik · yk

)(5.10)

für i = 1, ..., n. Da wir für die Eindeutigkeit bestimmt haben, dass die Elementeauf der Hauptdiagonlen von L jeweils gleich 1 sind, gilt:

yi = bi −i−1∑k=1

lik · yk (5.11)

5.1.1.3 Rückwärtseinsetzen

Abschließend muss noch das Gleichungssystem Ux = y gelöst werden. Das Vorge-hen ist analog zum Gauß-Verfahren und wird berechnet durch:

xi = 1uii

yi −n∑

k=i+1uik · xk

(5.12)

für i = n, ..., 1.Mit x erhalten wir dann die Lösung für das Gleichungssystem ax = b.

26

5 Lösen linearer Gleichungssysteme mit SkelCL

5.2 Umsetzung mit SkelCLWir sind schrittweise vorgegangen um das Lösen linearer Gleichungssysteme indem SkelCL Framework umzusetzen. Zunächst haben wir die einzelnen Schrit-te des Gauß-Verfahrens und der LU Dekomposition analysiert und eine Parallel-isierungsstrategie zu den einzelnen Schritten entwickelt. Anschließend haben wiruntersucht welche der Skeletons sich für die jeweiligen Algorithmen eignen könntenund haben die Algorithmen implementiert.

5.2.1 Vorwärtselimination/LU-ZerlegungDer Algorithmus der Vorwärtselimination des Gauß-Verfahrens und der Algorith-mus der LU-Zerlegung unterscheiden sich dadurch voneinander, das die Umfor-mungsschritte ebenfalls gespeichert werden.Beide Algorithmen arbeiten In-Place, sprich sie erzeugen keinen neuen Speicherin dem das Ergebnis steht sondern überschreiben den ursprünglichen Speicher. ImFalle der LU-Zerlegung werden sowohl die obere als auch die untere Dreiecksma-trix in einer Matrix gespeichert. Die Einsen auf der Hauptdiagonale von L werdendabei nicht gespeichert. Daraus ergibt sich folgende Anordnung im Speicher:

u1,1 u1,2 ... u1,n

l2,1 u2,2 ... u2,n... . . . . . . ...

ln,1 . . . ln,n−1 un,n

(5.13)

Im Falle des Gauß-Verfahrens wird die um b erweiterte Koeffizientenmatrix imSpeicher abgebildet. Siehe Formel 5.4

Die Vorwärtselimination ist wie bereits erwähnt ein iterativer Prozess und bestehtim Grunde daraus ein Vielfaches einer Zeile zu allen darunter liegenden Zeilen hin-zu zu addieren. Die einzelnen Elemente einer Zeile sind unabhängig voneinanderund können daher Parallel berechnet werden. Auch die einzelnen Zeilen sind inner-halb einer Iteration unabhängig voneinander. Daher können wie in Abbildung 5.1dargestellt, große Bereiche parallel berechnet werden. In der Ersten Iteration kannder blaue Bereich berechnet werden, in der zweiten Iteration der rote Bereich, etc.Allerdings wird ebenfalls dargestellt, dass der Bereich mit zunehmender Iterationimmer kleiner wird.

Umgesetzt wurde der Algorithmus mit einem map-Skeleton, da es das flexibelstealler Skeletons ist und als einziges der Skeletons ermöglicht die 1-D Position beiVektoren oder eine 2-D Position bei Matrizen für die innerhalb der Berechnung

27

5 Lösen linearer Gleichungssysteme mit SkelCL

Abbildung 5.1: Unabhängige Daten bei der Vorwärtselimination

bereitzustellen. Die aktuelle (x,y)-Position der Elemente innerhalb der Matrix istallerdings für die Berechnung notwendig da allgemein für die i-te Iteration undy > i, x ≥ i gilt:

uyx = uyx −uyi

uii

uix (5.14)

Das zu wählende Vielfache einer Zeile y ist also für die i-te Iteration:

− uyi

uii

(5.15)

und beruht damit auf einem Wert der zu transformierenden Zeile. Daher werdendiese Umformungsfaktoren im Vorfeld jeder Iteration separat, ebenfalls in einemmap-Skeleton berechnet und in einem Vektor zwischengespeichert.Im Falle der LU-Zerlegung werden die Umformungsfaktoren im Anschluss an Trans-formation zusätzlich in der Matrix gespeichert.

Listing 5.1: Vorwärtselimination mit SkelCL Skeletons.double gaussForwardMap(Matrix<float> &M){

size_t size = M.rowCount();//Definition des map-Skeletons fuer die VorwaertseliminationMap<void(IndexPoint)> forwardFunc(

"void func( IndexPoint position, float_matrix_t m, "\" __global float* lf,"\" int cols,int iteration)"\

28

5 Lösen linearer Gleichungssysteme mit SkelCL

"{ "\"position.x+=iteration; position.y+=iteration+1;"\"if(position.x>=cols) return;"\"float ret = get(m,position.y, position.x);"\"ret -= get(m, iteration, position.x)*lf[position.y];"\"set(m, position.y, position.x, ret);"\

"}");

//Definition des map-Skeletons fuer die Berechnung der UmformungsfaktorenMap<void(Index)> calcLineFactors(

"void func( Index i, float_matrix_t m, __global float* lf,int iteration)"\"{ "\

"lf[i] = get(m, i, iteration)/get(m, iteration, iteration);"\"}");

pvsutil::Timer timer; // ZeitmessungIndexVector index(size); // X-Position bereitstellenVector<float> lf(size); // Vektor fuer Umformungsfaktorenfor(int i=0; i< (int)size-1 ;i++){

calcLineFactors(index, M, lf, i); // berechnen der UmformungsfaktorenIndexMatrix pos({size-i, size+1-i}); // X,Y-Position bereitstellenforwardFunc(pos, M, lf, (int)size+1, i);// berechnen der i-ten Iteration

}M.dataOnDeviceModified(); // Daten wurden auf GPU geaendertM.copyDataToHost(); // Daten von der GPU kopierenreturn timer.stop();

}

In Listing 5.1 sehen Sie den Quelltext für die Vorwärtselimination des Gauß-Verfahrens. Es werden dort zunächst die zwei beschriebenen map-Skeletons mit-samt der Funktionsdefinition definiert. Anschließend wird in einer Schleife überdie einzelnen Iterationen iteriert und die Skeletons ausgeführt. Da der Algorith-mus In-Place arbeitet muss dem Framework abschließend noch mitgeteilt werden,dass die Daten auf der GPU durch den Algorithmus geändert wurden und diesewieder in den Arbeitsspeicher des Host kopiert werden sollen.Zu beachten ist ebenfalls der Anfang der Funktionsdefinition des forwardFunc-map-Skeletons, hier wird zunächst mit der Zeile "position.x+=iteration; positi-on.y+=iteration+1;" der Index nach unten-rechts verschoben. Das ist notwendig,da die verwendete IndexMatrix nur (x,y)-Koordinaten beginnend ab (0,0) liefert.Unser Algorithmus arbeitet aber wie Abbildung 5.1 zeigt auf einem Variablen Aus-schnitt der Matrix. Eine Weitere Besonderheit ist die Zeile ïf(position.x>=cols)return;". Bedingt durch die Architektur der GPUs werden in der Regel immer einefeste Anzahl von Threads auf der GPU ausgeführt, daher kann es vorkommen,dass obwohl wir durch die IndexMatrix die Größe des zu berechnenden Bereichsbestimmt haben, trotzdem ein größerer Bereich bearbeitet wird. Da die Daten Se-riell im Speicher liegen und die get() und set() Funktionen keinerlei Überprüfungenthalten ob die Grenzen eingehalten werden, muss diese Überprüfung manuellerfolgen. Diese Einschränkung trifft aber nur auf Algorithmen zu, die In-Place ar-beiten.

29

5 Lösen linearer Gleichungssysteme mit SkelCL

Abbildung 5.2: Unabhängige Daten bei dem Rückwärtseinsetzen

5.2.2 Vorwärtseinsetzen/RückwärtseinsetzenBetrachtet man die beiden Formeln 5.10 für das Vorwärtseinsetzen und 5.12 fürdas Rückwärtseinsetzen so erkennt man, dass sich beide Formeln sehr ähnlich sind.Sie unterschieden sich abgesehen von den Bezeichnungen dadurch, dass bei For-mel 5.10 die Summe über k = 1, ..., i − 1 gebildet wird und die Elemente in derReihenfolge i = 1, ..., n berechnet werden, wohingegen bei Formel 5.12 die Summeüber k = i + 1, ..., n gebildet wird und die Elemente in der Reihenfolge i = n, ..., 1berechnet werden. Das Problem ist also abgesehen von den Laufrichtungen iden-tisch. Daher wird an dieser Stelle zunächst nur das Rückwärtseinsetzen im Detailbeschrieben.

Auch diesen Algorithmus haben wir analysiert und auf unabhängige Daten hinuntersucht. Die Reihenfolge, in der die Elemente berechnet werden ist vorgegebenda der Wert von xn−1 von dem Wert von xn abhängt. Daraus folgt, dass es nichtmöglich ist die Berechnungen für x1, ..., xn naiv parallel ablaufen zu lassen.Es wäre zwar möglich die Berechnung der Summe zu parallelisieren, allerdings fielnach genauerer Analyse auf, dass der Algorithmus auch anders aufgebaut werdenkann. Statt der Reihe nach jeweils über ein einzelnes Element xi zu iterieren kannauch über die Elemente der Summe uikxk iteriert werden. Dadurch kann eine hö-here Parallelisierung erreicht werden als bei der parallelen Berechnung der Summe.

30

5 Lösen linearer Gleichungssysteme mit SkelCL

Betrachten wir zum besseren Verständnis das Beispiel aus Abbildung 5.2. Zu Be-ginn wird hier y8 = y8

u88errechnet. In der folgenden Iteration wird nun y7 = y7 −

u78y8, y6 = y6− u68y8, ..., y1 = y1− u18y8 gerechnet. In der Abbildung ist dies in inblau gekennzeichnet. Anschließend wird y7 = y7

u77gesetzt. In der folgenden Iterati-

on (rot gekennzeichnet) wird nun y6 = y6−u67y7, y5 = y5−u57y7, ..., y1 = y1−u17y7gesetzt.

Allgemein gilt folgende Prozedur:

Begonnen wird einmalig mit:yn = yn

unn

(5.16)

Danach wird für jede Iteration i = n− 1, ..., 1 in dieser Reihenfolge ausgeführt:

1.yj = yj − ujiyi,∀j ∈ {1, ..., n− 1} (5.17)

2.yi = yi

uii

(5.18)

Nachdem alle Iterationen durchlaufen sind, steht in dem Vektor y das gesuchteErgebnis. Der Algorithmus für das Vorwärtseinsetzen läuft analog ab, aber mitunterschiedlichen Laufrichtungen.Durch diese alternative Interpretation des Algorithmus gelangen wir zu dem glei-chen Endergebnis, erzielen aber eine höhere Parallelität. Wir haben den Algorith-mus so implementiert, dass die Parallelität noch weiter gesteigert werden kann,indem mehrere Lösungsvektoren gleichzeitig errechnet werden.

Auch dieser Algorithmus wurde mit zwei map-Skeletons implementiert wie Lis-ting 5.2 zeigt. Das erste map-Skeleton backwardFunc implementiert die Formel5.17 und das zweite map-Skeleton divideLine implementiert die Formel 5.18.Beide Skeletons werden in der Schleife n−1 mal aufgerufen und das backwardFuncwie es die Prozedur vorsieht einmal zusätzlich vor der Schleife. Auch hier tauchenwieder Abbruchbedingungen auf um die Grenzen des zu bearbeitenden Abschnittszu überprüfen.

Listing 5.2: Rückwärtseinsetzen mit SkelCL Skeletons.double gpuBackward(Matrix<float> &LU, Matrix<float> &Y, Matrix<float> &X){

size_t sizeR = LU.rowCount();size_t sizeC = Y.columnCount();Map<void(IndexPoint)> backwardFunc(

"void func( IndexPoint position, float_matrix_t lu,"\" float_matrix_t x, "\" size_t cols, size_t rows, int iteration )"\

31

5 Lösen linearer Gleichungssysteme mit SkelCL

"{ "\" if(position.x>=cols) return;"\" if(position.y>=(rows-iteration)) return;"\" float ret = get(x,position.y, position.x);"\" ret -= get(lu, position.y, rows-iteration) "\" * get(x,rows-iteration , position.x);"\" set(x, position.y, position.x, ret);"\"}");

Map<void(Index)> divideLine("void func( Index i, float_matrix_t lu, "\" float_matrix_t x, int line, size_t cols)"\" { "\" if(i>=cols) return;"\" float ret = get(x,line, i);"\" set(x, line, i, ret / get(lu,line, line) );"\" }");

pvsutil::Timer timer;IndexVector index(sizeC);

divideLine(index, LU, X, sizeR-1, sizeC);for(int i=1; i< (int)sizeR ;i++){

IndexMatrix pos({sizeR-i, sizeC });backwardFunc(pos, LU, X, sizeC, sizeR, i);divideLine(index, LU, X, sizeR-i-1, sizeC);

}

X.dataOnDeviceModified();X.copyDataToHost();return timer.stop();

}

5.2.3 NachiterationDas Gauß-Verfahren und die LU Dekomposition sind Verfahren, die auf sehr vielenAdditionen von Fließkommazahlen beruhen, daher können sich bereits kleine Run-dungsfehler sehr schnell aufschaukeln bei Matrizen großer Größen. Dadurch kön-nen erhebliche Abweichungen zwischen der CPU Implementierung und der GPUImplementierung entstehen. In unseren Tests konnten wir bei Berechnungen ineinfacher Genauigkeit und Matrizen der Größe 512x512 bereits Abweichungen von10−1 beobachten. Selbst bei doppelter Genauigkeit konnten wir bei unseren TestsAbweichungen von 10−11 beobachten. Die CPU rechnet intern mit 80 Bit und istdaher immer genauer als die GPU.

Eine Möglichkeit die Rechengenauigkeit zu erhöhen ist die Nachiteration. Dabeistartet man mit der berechneten Lösung für x als x0 und errechnet dann dasResidium:

rk = b− Axk (5.19)Anschließend berechnet man mit der bekannten LU-Zerlegung die Lösung für zk

Azk = rk (5.20)

und setzt abschließend:xk+1 = xk + zk. (5.21)

32

5 Lösen linearer Gleichungssysteme mit SkelCL

Dieser Prozess kann solange wiederholt werden bis man die gewünschte Genauig-keit erreicht hat. In der Praxis hat sich gezeigt das eine Iteration schon ausreichtum die gewünschte Genauigkeit zu erreichen.

Bei der Nachiteration zeigt das SkelCL Framework seine Stärken. Mit ein paarwenigen Zeilen Code wie, im Folgenden Listing dargestellt, lässt sich die Nachite-ration durchführen.Die ersten Zeilen definieren die Funktionen plus minus und eine Matrixmultipli-kation. Nach der Definition kann z.B. Formel 5.19 mit der einer einzelnen Zeileberechnet werden: Matrix<float> r = minus(B, matrixMult(M,X));.

Listing 5.3: Nachiteration mit SkelCL Skeletons.double gpuNachiteration(Matrix<float> &M, Matrix<float> &B, Matrix<float> &LU, Matrix<float> &X){

Zip<float(float, float)> matrixMultZip("float func(float x, float y){ return x*y; }");Reduce<float(float)> matrixMultReduce("float func(float x, float y){ return x+y; }");AllPairs<float(float, float)> matrixMult(matrixMultReduce, matrixMultZip);Zip<float(float, float)> minus("float func(float x, float y){ return x-y; }");Zip<float(float, float)> plus("float func(float x, float y){ return x+y; }");pvsutil::Timer timer;

Matrix<float> Y_temp({X.rowCount(), X.columnCount()});Matrix<float> X_temp({X.rowCount(), X.columnCount()});

Matrix<float> r = minus(B, matrixMult(M,X));printf("Residium max %f \n", getMaxMatrixBetrag(r));

//GPU forwardgpuForward(LU, r, Y_temp);

//GPU backwardgpuBackward(LU, Y_temp, X_temp);printf("X_temp max %f \n", getMaxMatrixBetrag(X_temp));

Matrix<float> xn = plus(X, X_temp);copy(xn.begin(), xn.end(), X.begin());X.dataOnHostModified();

Matrix<float> r2 = minus(B, matrixMult(M,X));printf("Residium max nach korrektur %f \n", getMaxMatrixBetrag(r2));

return timer.stop();}

5.3 ErgebnisseZiel war es das Lösen linearer Gleichungssysteme durch den Einsatz des SkelCLFrameworks und somit von parallelen Berechnungen auf der GPU zu beschleuni-gen. Dazu haben wir den Algorithmus ebenfalls auf der CPU implementiert unddie Laufzeiten miteinander verglichen. Die Ergebnisse sind in der folgenden Ta-belle aufgelistet. Wie Sie sehen können wurde mit der Intel HD 4000 GPU eineBeschleunigung von bis zu 23x erreicht werden und mit der NVIDIA GeForce GTX

33

5 Lösen linearer Gleichungssysteme mit SkelCL

TITAN GPU sogar eine Beschleunigung von bis zu 510x. Jeweils im Verhältnis zurLaufzeit der Single-Core CPU Implementierung auf einem 2,5 GHz Intel Core i5.Man sieht aber auch, dass sich der Einsatz von GPUs besonders lohnt, wenn großeDatenmengen verarbeitet werden sollen. Bei Geringer Problemgröße der Matrixbzw. Anzahl der zu berechnenden Lösungen war die CPU Implementation bis zuzweimal so schnell wie die GPU Implementierung.

Matrix CPU Intel HD 4000 GeForce GTX TITAN(2,5 GHz Intel Core i5) (einfache Genauigkeit) (doppelte Genauigkeit)

(nxn)(nxn) Laufzeit Laufzeit SpeedUp Laufzeit SpeedUp64 26 ms 50 ms 0,25x128 222 ms 104 ms 2,13x256 1.917 ms 229 ms 8,37x512 15.905 ms 706 ms 22,52x 144 ms 110,45x1024 143.608 ms 6.137 ms 23,40x 576 ms 249,32x2048 1.398.900 ms 60.479 ms 23,13x 2.741 ms 510,36x

34

6 Fazit und AusblickDas Masterprojekt und die Zusammenarbeit mit den beiden Masterarbeiten boteinen interessanten Einstieg in die Welt der parallelen Algorithmen und GPG-PU Computing. Nach einer gewissen Lernphase konnten erste Problemstellungenelegant und effizient gelöst werden. Für den maximalen Nutzen der GPUs sind je-doch viele feingranulierte Optimierungen nötig, die zunächst viel Arbeit und Zeitkosten. Dies mag jedoch auch daran liegen, dass uns die Programmierung vonCPUs vertrauter ist und muss keinen prinzipiellen Nachteil der GPGPU Program-mierung darstellen. Systeme, bei denen mögliche Parallelisierungen nicht direktklar ersichtlich sind erfordern komplexe Analysen und kreative Lösungsstrategi-en. Selbst High-Level Abstraktionen, wie das vorgestellte SkelCL, verlieren beinicht trivialen Herausforderungen schnell an Handlichkeit. Es kann jedoch festge-halten werden, dass bei sehr großen Datenmengen ein Einsatz von GPUs massiveGeschwindigkeitsvorteile bringen kann. Um diese Vorteile allerdings optimal aus-nutzen zu können müssen zunächst die Algorithmen angepasst und zum großenTeil komplett neu konzipiert werden. Viele etablierte Optimierungen von Algo-rithmen sind schlicht nicht für eine Parallelisierung geeignet. Es muss daher oftnach performanten Alternativen gesucht werden.

Die Erkenntnis des Projekts ist, dass in vielen Bereichen der Informatik paral-lele Algorithmen und GPGPU Computing eine sehr große Rolle spielen und auchin Zukunft eine große Rolle spielen werden. Die Anwendungen werden immer kom-plexer und rechenintensiver. Dies gilt insbesondere auch für die Forschung. Es isthäufig nur möglich den aktuellen und zukünftigen Anforderungen durch die Nut-zung paralleler Rechner und Rechencluster gerecht zu werden. Produktive Ent-wickler sollten aber bereits über gewisse Erfahrungen in diesem Feld verfügen, dasonst der Kosten-Nutzen Faktor sehr schnell ins negative ausschlagen kann. DiesesForschungsfeld ist ein noch vergleichbar junges, weshalb noch mit einigen Fort-schritten zu rechnen ist. Dies bezieht sich sowohl auf die Fähigkeiten der GPUsals auch auf die Frameworks und APIs für die Programmierer.

35

Literaturverzeichnis[1] Cuda c programming guide, Feb. 2014. http://docs.nvidia.com/cuda/

cuda-c-programming-guide/index.html.

[2] Was ist gpu computing, Jan. 2014. http://www.nvidia.de/object/gpu-computing-de.html.

[3] R. Böing. Implementation eines cuda basierten kalman-filters zur spurrekon-struktion des atlas-detektors am lhc, Feb. 2014. http://www.lab4inf.fh-muenster.de/lab4inf/docs/thesis/MA_Rene_Boeing.pdf.

[4] R. Frühwirth. Application of Kalman filtering to track and vertex fitting. Nucl.Instrum. Methods Phys. Res., A, 262(HEPHY-PUB-503):444. 19 p, Jun 1987.

[5] R. E. Kalman. A new approach to linear filtering and prediction problems.Transactions of the ASME-Journal of Basic Engineering, 82(Series D):35–45,1960.

[6] J. Pequenao and P. Schaffner. An computer generated image representing howATLAS detects particles. http://cds.cern.ch/record/1505342?ln=de, Jan 2013.

[7] M. Steuwer and S. Gorlatch. High-Level Programming for Medical Imaging onMulti-GPU Systems using the SkelCL Library. 2International Conference onComputational Science, ICCS 2013, (DOI 10.1016/j.procs.2013.05.239), Jun2013.

[8] M. Steuwer, P. Kegel, and S. Gorlatch. SkelCL – A Portable Skeleton Libra-ry for High-Level GPU Programming. 2011 IEEE International Parallel andDistributed Processing Symposium, (DOI 10.1109/IPDPS.2011.269), Jun 2011.

[9] G. Welch and G. Bishop. An introduction to the kalman filter. Technicalreport, Chapel Hill, NC, USA, 2006. http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf.

36