Parallelle abstracties voor het programmeren van GPU's in Julia

Parallelle abstracties voor het programmeren van GPU's in Julia Pieter Verstraete

Promotor: prof. dr. ir. Bjorn De Sutter Begeleider: Tim Besard Masterproef ingediend tot het behalen van de academische graad van Master of Science in de ingenieurswetenschappen: computerwetenschappen

Vakgroep Elektronica en Informatiesystemen Voorzitter: prof. dr. ir. Jan Van Campenhout Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2013-2014

Voorwoord GPU’s zijn de laatste jaren uitgegroeid tot zeer krachtige hardwarecomponenten, gespecialiseerd in parallelle verwerking van grote hoeveelheden data. In veel wetenschappelijk onderzoek levert het gebruik van een GPU dan ook een grote prestatiewinst op bij het uitvoeren van zware rekentaken. Het programmeren van een GPU is echter complex en een algoritme voor een GPU implementeren vereist technische kennis en veel ontwikkelingstijd van de programmeur. In veel onderzoeksgroepen beschikt men niet over deze kennis en wil men zeer snel algoritmes kunnen implementeren, testen en evalueren. In deze masterproef heb ik onderzoek kunnen verrichten in dit boeiend en belangrijk onderzoeksgebied. Bovendien combineert dit onderzoek een aantal van mijn grote interesses binnen het vakgebied computerwetenschappen, namelijk GPU’s, compilers en productiviteit van softwareontwikkeling. Ik heb dit project dan ook met veel plezier gedaan. Het is mij een aangename plicht alle mensen te bedanken die een bijdrage hebben geleverd tot het welslagen van deze masterproef. In de eerste plaats dank ik mijn promotor, prof. dr. ir. Bjorn De Sutter, voor deze opportuniteit en zijn visie op dit onderwerp. Verder dank ik mijn dagelijkse begeleider, Tim Besard, voor zijn voortdurende steun. Deze masterproef was voor hem een hoofdbekommernis en hij heeft vele uren vrijgemaakt om mij met allerhande problemen te helpen. Tenslotte dank ik ook ouders, familie en vrienden voor de steun en aangename momenten de voorbije jaren.

Pieter Verstraete, juni 2014

i

Toelating tot bruikleen De auteur(s) geeft(geven) de toelating deze masterproef voor consultatie beschikbaar te stellen en delen van de masterproef te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze masterproef.

Pieter Verstraete, juni 2014

ii

iii

Parallelle abstracties voor het programmeren van GPU’s in Julia door Pieter Verstraete Promotor: prof. dr. ir. Bjorn De Sutter Begeleider: Tim Besard Masterproef ingediend tot het behalen van de academische graad van Master of Science in de ingenieurswetenschappen: computerwetenschappen Vakgroep Elektronica en Informatiesystemen Voorzitter: prof. dr. ir. Jan Van Campenhout Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2013-2014

Samenvatting Deze masterproef onderzoekt de mogelijkheid om de rekenkracht van een GPU beschikbaar te maken in de hoogniveau programmeertaal Julia. Een belangrijke vereiste hierbij is dat de hoge productiviteit van de programmeur die Julia biedt, gewaarborgd blijft. Om een GPU te programmeren is immers een zeer gespecialiseerde en gedetailleerde technische kennis vereist, en dit maakt het werk moeilijk en tijdrovend. In de moderne wetenschap is vaak enorm veel rekenkracht nodig, maar juist in die middens ontbreekt vaak deze technische kennis om de kracht van GPU’s te benutten. Het doel van deze masterproef is de GPU toegankelijk te maken in Julia op een manier die de details van de GPU afschermt van de programmeur. In een eerste stap definiëren en implementeren we abstracties in Julia die in parallel op de GPU kunnen uitgevoerd worden. Vervolgens hebben we de Julia compiler aangepast om deze abstracties te vertalen naar GPU instructies. Ook is een infrastructuur ge¨ımplementeerd die de GPU beheert op een manier die voor de programmeur transparant is. Tenslotte worden de abstracties en ontwikkelde infrastructuur geëvalueerd aan de hand van een concrete toepassing, namelijk de tracetransformatie.

iv

Trefwoorden GPU, Julia, parallelle abstracties, LLVM, PTX

Parallel abstractions for programming GPUs in Julia Pieter Verstraete Supervisor(s): Bjorn De Sutter, Tim Besard Abstract— This master’s thesis explores the possibility to provide access to the computing power of a GPU from the high-level programming language Julia. An important requirement here is to keep the programmer’s productivity at the same high level as if he would use Julia without a GPU. Indeed, very specialized and detailed technical knowledge is needed in order to program a GPU, making it complex and time-consuming. In many modern scientific domains quite a lot of brute computing power is required, but often these domains lack the technical expertise to use GPUs in an efficient manner. The purpose of this thesis is to provide access to a GPU from Julia in a way that shields the GPU details from the programmer. In a first step we define and implement in Julia abstractions that can be executed in parallel on the GPU. Next we adapt the Julia compiler such that it can translate these abstractions to GPU code. The resulting compiler infrastructure manages the GPU in a way that is transparent to the programmer. Finally we evaluate the abstractions and compiler infrastructure in the context of a concrete application, namely the trace transform. Keywords—GPU, Julia, parallel abstractions, LLVM, PTX

I. I NTRODUCTION GPUs are very powerful computers specialized in the execution of parallel mathematical calculations. However, tapping this computing power is in many scientific domains hampered by the lack of sufficient technical knowledge needed for programming a GPU. Indeed, most software is developed in high-level programming languages because they offer the necessary abstractions and thereby allow the programmer to work efficiently. However, high-level languages usually cannot provide sufficient execution performance. One high-level programming language of interest in this thesis is Julia. It is a dynamic programming language that is used extensively in mathematical and scientific contexts. It enables very high programmer productivity but it currently lacks an infrastructure to access a GPU. The challenge addressed in this work is to provide access to a GPU from the language Julia in such a way that shields the GPU details from the programmer. The objective is to combine high programmer productivity with efficient program execution on the GPU. Our approach is to define and implement a number of abstractions in Julia. Doing the implementation in Julia allows analysis and optimization based on runtime information. Furthermore, the implementation framework takes care of most of the administrative tasks needed in using GPUs, such as loading GPU code, allocation and transport of memory, etc. The resulting compiler framework is evaluated using a specific case, namely the trace transform. This is an application that is inherently highly parallel and lends itself very well to the use of a GPU. Both programmer productivity and execution time are evaluated.

II. R ELATED WORK GPUs have been used for a long time in the context of computer graphics and imaging. In particular, computer games make heavy use of GPUs. More recently, a more general broader scientific community began to show interest in GPUs because of its enormous parallel computing power. Therefore, a lot of research was done to find ways to tap the power of a GPU in a way that was both effective and efficient. The first step in this direction was CUDA, along with CUDA C, which provide a more general access to the graphical pipeline of the GPU. But in spite of the very good execution performance of the CUDA code, its use is very complex, time-consuming, and error-prone. In this paper we have chosen to build abstractions in Julia [4] that provide access to the GPU in a productive but efficient manner. Abstractions have been shown to improve the productivity of software development [1], [2]. Our strategy is similar to [3], but differs in two major aspects. First of all, the code is compiled dynamically such that runtime information can be used to optimize the abstractions. And secondly, our compiler framework also takes care of the administrative tasks required for this access to the GPU. In order to evaluate the results of our work we used a concrete application: the trace transform [5]. This mathematical function extracts certain properties from an image and is used in numerous domains such as object recognition, facial recognition, and so on. Although this transform is very robust and powerful, it is not often applied because of the high computational complexity. Since the trace transform is inherently very parallel, implementing it on a GPU gives a tremendous execution time improvement. Hence, this transform is a good way to evaluate the results of our work. III. C OMPILATION OF J ULIA TO PTX As a first step in our approach, we create a compiler infrastructure that meets two requirements: • A programmer can write GPU kernels in the Julia language. This removes the need for a two-phase development. • The administrative tasks of using a GPU are taken care of by the compiler infrastructure, which increases programmer productivity. A. API in Julia Writing GPU kernels in Julia is similar to writing a regular function, apart from one difference. A GPU reads from and writes to memory and no actual values are explicitly returned. Hence, a Julia kernel function must always explicitly return a nothing value and declare output variables as function parameters.

Furthermore, writing GPU kernels requires access to some extra functionality, which has been added to Julia as an API: • Shared memory can be defined using shmem = cuSharedMem() and accessed using setCuSharedMem(shmem, index, value) and getCuSharedMem(shmem, index). • The intrinsics for the thread and block id, for the block and grid dimension, and for thread synchronization are accessible via simple function calls. We defined a macro which performs the administrative tasks of using a GPU. Firstly, this macro will allocate memory on the GPU for both input and output parameters and transport input values to the GPU. Then, function compilation is triggered and the resulting PTX code is loaded in the GPU. After this, the GPU is signaled to execute the kernel. When execution is finished, the output values are copied back to the CPU memory and the allocated memory is freed. B. Modifications to the Julia compiler The Julia compiler infrastructure is based on LLVM. A back end for compiling to PTX instructions is available, namely NVPTX, but it does not fully support the complete LLVM IR. Hence, in our work we re-used the Julia front end but we adapted the code generation step in order to address the following restrictions. • Dedicated modules: The code that is targeted for a GPU must be wrapped in a dedicated LLVM module, which is then compiled into PTX code. Hence the programmer must indicate which functions are to be implemented as GPU kernels. • Host vs. device functions: A host function is launched by the CPU to access the GPU kernel; a device function can only be called from the GPU code. The NVPTX back end requires that the LLVM IR functions which are meant to be host functions have an extra annotation. • NVPTX intrinsics: CUDA C offers a number of built-in variables and functions in the form of intrinsics that can be used by the LLVM IR. These intrinsics must also be offered in Julia. • Alternative numbering: Array indexes in Julia start from 1, while in C they start from 0. Also indexes for threads and blocks in Julia start from 1. Hence, a transparent mapping back-and-forth must be foreseen by the compiler infrastructure. • Memory management: CUDA supports different types of memory, most notably global and shared memory, each with their own address space. Julia, however, only supports the default address space. The compiler infrastructure that we developed must therefore provide a mapping from the Julia memory management to the CUDA memory concepts. Furthermore, the compiler must be adapted to allow specifying which function arguments are input and which are output. Finally, array bound checking is automatically enforced in Julia, while CUDA offers no such notion. • Link with libdevice: NVIDIA offers libraries of mathematical functions, written in LLVM bitcode, but these are optimized for a specific GPU architecture. Hence, the correct libdevice version must be selected based on the target harware. This selection and linking process is done by the

CUDA compiler and must now also be done in the Julia compiler infrastructure. IV. A BSTRACTIONS When choosing which abstractions would be provided to the programmer, we took into account two things. Firstly, the abstractions were to be executed on a GPU and thus had to exhibit some parallel properties. Secondly, as the trace transform is used for our evaluation, the abstractions had to be useful for this calculation, but of course also still generically useful. Based on these guidelines and [7], the following well known abstractions were chosen: map, scan and reduce. We chose to make use of the infrastructure described above and implement these abstractions in Julia for three reasons. First of all, the fact that Julia code is JIT compiled allows analyzing and optimizing the abstractions code based on runtime information, which is not possible when statically compiling CUDA C code. Furthermore, the definition and implementation of the abstractions is decoupled from the CUDA framework. When providing support for different infrastructures, like OpenCL, the existing implementation of the abstractions can be reused. Finally, Julia as a high-level language allows the programmer to implement these and future abstractions in a productive manner. These abstractions can now be used to implement the T1 functional, which can be seen in Equation (1). T 1(column) =

# rowsX in column i=m

column[i] ∗ (i − m)

(1)

Here, m is the index of the weighted median of the column. Calculation of this weighted median happens in two steps. First the column is scanned in order to find the total and intermediary sums of the elements. Then the median index is the value i for which scanned column[i − 1] < total 2 . This constraint can be and scanned column[i] ≥ total 2 checked for all positions independently, hence a map is used. The calculations of Equation (1) are also implemented using two abstractions. First a map is applied which multiplies all column elements with (i − m). The results of this map are then summed using a reduce. V. R ESULTS AND ANALYSIS Based on the infrastructure and abstractions described above, this section evaluates the improvement of the trace transform implemented in Julia by using a GPU, but without compromising the programmer’s productivity. The evaluation is done in multiple steps, starting from the best execution performance but lowest productivity. Successive steps increase the performance with as little impact as possible on productivity. A. Description of the experiments Each experiment implements on the GPU a rotation (from 1 to 359 degrees in steps of 1 degree) of the blackand-white image shown in Figure 1, followed by a T1 functional. We used an Intel i7-3770K (Quad core @ 3.5GHz) CPU and an NVIDIA GeForce GTX Titan GPU. The OS is

Kernel rotate prescan find weighted median t1

Relative execution time 1.52 0.98 1.74 1.08

TABLE II R ELATIVE EXECUTION TIME OF THE J ULIA KERNELS WITH RESPECT TO THE CUDA KERNELS

Execution time (s)

Fig. 1. Reference image used for evaluation

10−1

Abstractions Julia → PTX Julia & CUDA C++ & CUDA

10−2 64x64

128x128 256x256 512x512 Image resolution (pixels)

Fig. 2. Average execution time of the trace transform using the T1 functional for different image resolutions

Debian (jessie) based on Linux kernel 3.13-1-amd64 x86 64. The version of Julia used is 0.3.0-prerelease+3033 (2014-05-17 15:19 UTC). The version of CUDA is 5.5 and the LLVM is version 3.3. B. Performance analysis Table I summarizes the results. The following implementation methods were evaluated: B.1 Manual integration of CUDA code in Julia This approach makes use of the package CUDA.jl [6] and the same CUDA code as in the reference implementation in C++ and CUDA. We see little difference in performance. The only point worth making is the startup of the kernel: CUDA.jl is on average 60% slower than the native CUDA API. The reason is that all arguments for the kernel must be cast to a void** type, which requires encapsulation in Julia and is therefore slower. But the average impact of this on the overall execution time decreases as the size of the image increases; the CPU is able to encapsulate the variables and invoke kernels faster than the GPU can process them, which can be seen in Figure 2. B.2 Compilation of Julia to PTX The infrastructure for compiling Julia code to PTX was described in Section III. The results are shown in Table I. The reasons for the lower performance are twofold: • Execution time of the kernels themselves, which can be seen in Table II: Julia automatically chooses a 64-bit float-

ing point division while the baseline CUDA implementation makes use of 32-bit divisions. The programmer could force a 32-bit division but this goes against the principle of shielding details of the underlying hardware component. • Execution time of the CPU vs GPU: The @cuda macro in Julia takes care of quite some administrative tasks (memory allocation, transfer of data, retrieval of cached code) which impose overhead. This overhead is now substantial enough such that the GPU is often idle waiting for new code from the CPU. Referring to Figure 2 we see that if the resolution is low (64x64 or 128x128) then the CPU cannot keep up with the GPU and the performance is constrained by the (slow) speed of the CPU. B.3 Use of abstractions Table I shows the performance using abstractions. Similar to the previous case, also this time the @cuda macro is used, imposing a lower-bound on the performance for small images. The use of abstractions, however, requires four abstractions for calculating the T1 functional, while the previous case only uses three kernels. This extra kernel explains the 25% performance penalty. Obviously, one could combine two abstractions into a single one as they are executed one after the other anyway. But this means that one has to make available many more abstractions (combinations of multiple lower abstractions), which means more complexity and therefore lower productivity. Alternatively, the compiler could be made more intelligent and make use of global analysis to detect that the two abstractions could be combined into one. Using the JIT compilation approach of Julia, the generated code could then be adapted to take care of this. C. Productivity analysis How does one measure productivity? One possible indicator is the number of lines of code required for a given algorithm. If a given algorithm can be implemented with fewer lines of code then these lines are semantically richer and the programmer is more productive. Table III shows the number of lines of code needed for implementing the T1 functional. A two-step implementation in Julia and CUDA is already an improvement over C++ and CUDA because of the highlevel nature of Julia itself. Using the Julia compiler infrastructure developed in this master’s thesis one can further improve the productivity because now no CUDA code is needed. Still, implementing everything completely in Julia without a GPU is even more productive because all code is now done in one and the same high-level language. If the programmer makes use of abstractions, as described earlier,

Implementation environment Julia C++ & CUDA Julia & CUDA Julia → PTX Abstractions

Average execution time (s) 0.6668 0.0201 0.0203 0.0248 0.0302

Average relative speedup 1.0 33.17 32.85 26.89 22.08

TABLE I AVERAGE EXECUTION TIME OF THE TRACE TRANSFORM USING THE T1 FUNCTIONAL AND AVERAGE RELATIVE SPEEDUP WITH RESPECT TO THE J ULIA IMPLEMENTATION

Implementation environment Julia C++ & CUDA Julia & CUDA Julia → PTX Abstractions

T1 (SLOC) 19 147 172 96 20

TABLE III S OURCE LINES OF CODE (SLOC) NEEDED FOR IMPLEMENTING THE TRACE TRANSFORM , R ADON AND T1 FUNCTIONAL

then he does not have to worry about the technical aspects of the GPU but rather he can focus completely on the conceptual nature of the problem. As a result, the productivity is comparable to Julia without a GPU. D. Future work The results presented above show that a combination of productivity and performance is achievable, but still two observations suggest that there is room for improvement: • Data transport between CPU and GPU: The transport of data between CPU and GPU forms an important bottleneck. Hence, it is important to avoid unnecessary transfers. The @cuda macro always transfers all input data from CPU to GPU and then all output data from GPU to CPU. However, if the data needs to be passed on from one kernel to another then transfer to the CPU is not necessary. For example, the output of the rotation is used only by the T functional, which is also executed on the GPU. • Function executed in the abstraction: The programmer has to specify via an extra input parameter which function is to be called by an abstraction. But this seriously limits the programmer’s flexibility. A better approach would be to support the programmer to pass an anonymous function (e.g. (a,b)->a+b), which can then be in-lined in the abstraction and would incur little overhead. VI. C ONCLUSION In this thesis we have developed a compiler framework for Julia that offers the programmer a framework to access the computing power of a GPU while still maintaining the high programming productivity that is characteristic for high-level programming languages such as Julia. The essence of the approach is to define and implement abstractions that can be executed on the GPU in a manner that is transparent to the programmer. The implementation of the abstractions is done in Julia such that it is possible to adapt the generated code based

on information gathered at execution-time. Since the compilation of Julia to GPU instructions is not supported in the current Julia compiler, we extended it with an infrastructure to do so. Furthermore, a number of administrative tasks required for the use of the GPU have been included in the compiler framework. This provides again more shielding of the low-level details from the programmer. The results in Section V show that the performance of the trace-transform is about 50% slower when using our abstractions and infrastructure as compared to manually programming the GPU. Nevertheless, the productivity of the programmer using this framework has improved several times as compared to a manual implementation. In addition to that, the execution performance is 22 times better with this framework than with the native Julia that does not make use of a GPU. These results show that the chosen solution is rather promising. The implementation is still not yet complete but our framework in the state it is now can be seen as a prototype that explores the possibilities. This prototype can be used for further research to find other abstractions or find ways to improve the implementation of the abstractions, e.g. by making use of global optimization. R EFERENCES [1] Kurt Keutzer and Tim Mattson A Design Pattern Language for Engineering (Parallel) Software Intel Technology Journal, vol. 13, no. 4, 2009 [2] Asanovic, Krste and Bodik, Rastislav and Demmel, James and Keaveny, Tony and Keutzer, Kurt and Kubiatowicz, John and Morgan, Nelson and Patterson, David and Sen, Koushik and Wawrzynek, John and others A view of the parallel computing landscape Communications of the ACM, vol. 52, no. 10, pp. 56-67, 2009 [3] Eric Holk, Milinda Pathirage, Arun Chauhan, Andrew Lumsdaine and Nicholas Matsakis GPU Programming in Rust: Implementing HighLevel Abstractions in a Systems-Level Language IEEE 27th International Parallel and Distributed Processing Symposium Workshops & PhD Forum (IPDPSW), pp. 315-324, 2013 [4] Jeff Bezanson, Stefan Karpinski, Viral B. Shah and Alan Edelman Julia: A Fast Dynamic Language for Technical Computing CoRR, vol. abs/1209.5145 [5] Alexander Kadyrov and Maria Petrou, The Trace Transform and Its Applications, IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 23, no. 8, pp. 811-828, August 2001 [6] Lin Dahua, CUDA.jl, https://github.com/JuliaGPU/CUDA.jl [7] Michael McCool, James Reinders and Arch Robison Structured parallel programming: patterns for efficient computation

Inhoudsopgave Voorwoord

i

Toelating tot bruikleen

ii

Overzicht

iii

Extended abstract

iv

Inhoudsopgave

ix

Gebruikte afkortingen

xii

1 Inleiding

1

1.1

Probleemstelling

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Doelstelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

Structuur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

2 Gerelateerd werk 2.1

4

Grafische processor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.1.1

Het programmeren van een GPU met CUDA . . . . . . . . . . . .

5

2.1.2

Compilatieproces en interactie met de GPU . . . . . . . . . . . .

8

2.1.3

GPU architectuur . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

Productiviteit versus prestatie . . . . . . . . . . . . . . . . . . . . . . . .

11

2.3

Julia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.3.1

De taal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.3.2

Werking van de Julia compiler . . . . . . . . . . . . . . . . . . . .

14

Tracetransformatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.4

ix

x

INHOUDSOPGAVE

3 Compilatie van Julia naar PTX

18

3.1

API in Julia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

3.2

Beperkingen van de Julia compiler

. . . . . . . . . . . . . . . . . . . . .

22

3.2.1

Gescheiden codebeheer . . . . . . . . . . . . . . . . . . . . . . . .

23

3.2.2

Onderscheid tussen host en device functies . . . . . . . . . . . . .

23

3.2.3

NVPTX Intrinsics

. . . . . . . . . . . . . . . . . . . . . . . . . .

23

3.2.4

Alternatieve telling . . . . . . . . . . . . . . . . . . . . . . . . . .

23

3.2.5

Geheugenbeheer . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.2.6

Linken met libdevice . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.2.7

Een in Julia geschreven kernel opstarten . . . . . . . . . . . . . .

25

Aanpassingen aan de Julia compiler . . . . . . . . . . . . . . . . . . . . .

25

3.3.1

Gescheiden codebeheer . . . . . . . . . . . . . . . . . . . . . . . .

25

3.3.2

Onderscheid tussen host en device functies . . . . . . . . . . . . .

26

3.3.3

NVPTX Intrinsics

. . . . . . . . . . . . . . . . . . . . . . . . . .

26

3.3.4

Alternatieve telling . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3.5

Geheugenbeheer . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.3.6

Linken met libdevice . . . . . . . . . . . . . . . . . . . . . . . . .

29

3.3.7

Een in Julia geschreven kernel opstarten . . . . . . . . . . . . . .

30

3.3

4 Abstracties

32

4.1

Bestaande implementaties van de tracetransformatie . . . . . . . . . . . .

32

4.2

Keuze van abstracties . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.2.1

Map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

4.2.2

Reduce . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

4.2.3

Scan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.3

Implementatie van de abstracties . . . . . . . . . . . . . . . . . . . . . .

39

4.4

Gebruik van de abstracties in de tracetransformatie . . . . . . . . . . . .

42

5 Resultaten en analyse

44

5.1

Informatie over de opstelling . . . . . . . . . . . . . . . . . . . . . . . . .

44

5.2

Prestatieanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

5.2.1

Manueel gebruik van CUDA in Julia . . . . . . . . . . . . . . . .

47

5.2.2

Compilatie van Julia naar PTX . . . . . . . . . . . . . . . . . . .

48

5.2.3

Gebruik van abstracties . . . . . . . . . . . . . . . . . . . . . . .

50

5.3

Productiviteitsanalyse . . . . . . . . . . . . . . . . . . . . . . . . . . . .

51

5.4

Toekomstig werk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

INHOUDSOPGAVE

xi

6 Conclusie

55

Lijst van figuren

60

Lijst van tabellen

61

Gebruikte afkortingen ALU

Arithmetic Logic Unit

API

Application Programming Interface

AST

Abstract Syntax Tree

CPU

Central Processing Unit

CUDA

Compute Unified Device Architecture

FPU

Floating Point Unit

GPGPU

General Purpose GPU

GPU

Graphical Processing Unit

IR

Intermediaire Representatie

ISA

Instructieset Architectuur

JIT

Just In Time

LLVM

Low Level Virtual Machine

MEX

MATLAB Executable

PTX

Parallel Thread Execution

SFU

Special Function Unit

SIMT

Single Instruction Multiple Thread

SLOC

Source Lines Of Code

SM

Streaming Multiprocessor

xii

Hoofdstuk 1 Inleiding 1.1

Probleemstelling

Dankzij de grote populariteit van computerspellen en andere grafische toepassingen, zijn GPU’s uitgegroeid tot zeer krachtige hardware componenten, gespecialiseerd in het uitvoeren van sterk parallelle rekentaken. Dit type rekentaak komt echter niet alleen voor in grafische context. Het gebruik van een GPU voor algemeen rekenwerk, ook wel General Purpose GPU (GPGPU) gebruik, kan in zeer veel wetenschappelijk onderzoek een grote snelheidswinst opleveren. Voorbeelden gaan van kankeronderzoek[1] tot het versnellen van neurale netwerken[2] en weermodellen[3]. Een probleem bij dit algemene gebruik van een GPU is dat het programmeren van een GPU gebeurt aan de hand van een laagniveau taal. Net als bij algemene laagniveau talen wordt de programmeur geconfronteerd met architecturale en technische details. In veel situaties ontbreekt de nodige kennis of tijd om met dergelijke laagniveau talen te werken. Daarom wordt tegenwoordig steeds vaker overgeschakeld op hoogniveau talen, zoals MATLAB en Python. Deze talen bieden allerhande abstracties aan om de productiviteit van de programmeur te verhogen. Daartegenover staat dat de code geschreven in een hoogniveau taal vaak minder efficiënt is dan wat met een laagniveau taal mogelijk is, wat resulteert in een significant prestatieverschil. Dit impliceert dus dat in veel situaties een keuze gemaakt moet worden tussen prestatie van de code en productiviteit van de programmeur. De uitdaging is dus een omgeving te ontwerpen die vermijdt dat deze keuze gemaakt moet worden. Een mogelijke oplossing is Julia[4], een hoogniveau, dynamische programmeertaal die voornamelijk gebruikt wordt in wiskundige en wetenschappelijke context. Julia is ontwikkeld als antwoord op deze vraag en probeert de productiviteit te bieden van een hoogniveau taal, zoals MATLAB, terwijl toch de prestaties van een laagniveau taal, 1

2

HOOFDSTUK 1. INLEIDING

zoals C, benaderd of zelfs geëvenaard worden. Hoewel Julia in de praktijk vaak aan deze verwachtingen voldoet, geldt dit niet voor GPGPU ontwikkeling. Aangezien hiervoor in Julia geen ondersteuning is ingebouwd, wordt de programmeur bij het programmeren van een GPU opnieuw geconfronteerd met technische details, wat de productiviteit van Julia tenietdoet. Het doel van deze masterproef is te onderzoeken op welke manier de prestaties van een GPU beschikbaar gemaakt kunnen worden in Julia zonder daarbij de productiviteit van de programmeur te schaden.

1.2

Doelstelling

Een veelgebruikte methode om de productiviteit van de programmeur te verhogen is gebruik te maken van abstracties. Het doel van een abstractie is namelijk het aanbieden van een bepaald concept of functionaliteit, terwijl de technische implementatie ervan verborgen blijft. Dit laat de programmeur toe om zich op een hoger abstractieniveau bezig te houden met de concepten van algoritmes, zonder zich te hoeven buigen over de technische implementatiedetails. Voor dit onderzoek zullen dus een aantal abstracties gedefinieerd worden die het programmeren van een GPU op productievere wijze laat verlopen. De implementatie van deze abstracties zal gebeuren in de dynamische, hoogniveau programmeertaal Julia. Dit heeft meerdere voordelen, waaronder de mogelijkheid tot het uitvoeren van analyses en optimalisaties op basis van runtime informatie, wat bij statische compilatie niet kan. Dit vereist dan wel dat Julia code gecompileerd moet kunnen worden naar GPU instructies. In deze masterproef zal dus enerzijds de Julia compiler uitgebreid worden om deze compilatie mogelijk te maken. Daarnaast zal ook een pakket ge¨ımplementeerd worden dat een aantal administratieve taken in verband met het gebruik van een GPU op zich neemt, zoals het inladen van de GPU code, allocatie en transport van geheugen, etc. Beide aspecten hebben tot doel de productiviteit van de programmeur te verhogen zonder aan performantie in te boeten. Voor de evaluatie van de behaalde resultaten zal gebruik gemaakt worden van een specifieke toepassing, namelijk die van de tracetransformatie. Dit is een algoritme dat probeert getransformeerde versies van dezelfde figuur te herkennen. Dit algoritme bevat sterk parallelle taken en leent zich dus zeer goed tot het gebruiken van een GPU. Aangezien deze masterproef gaat over het verenigen van prestatie en productiviteit zullen beiden geëvalueerd worden.

1.3. STRUCTUUR

1.3

3

Structuur

Na in dit hoofdstuk de probleemstelling en de oplossing vermeld te hebben, wordt in Hoofdstuk 2 dieper ingegaan op de nodige technische achtergrond. In Hoofdstuk 3 wordt de infrastructuur voor het compileren van Julia code naar GPU instructies en het uitvoeren van administratieve taken, uit de doeken gedaan. Vervolgens wordt de definitie en implementatie van de gekozen abstracties besproken in Hoofdstuk 4. De analyse en evaluatie van de bekomen resultaten gebeurt in Hoofdstuk 5. Tenslotte bevat Hoofdstuk 6 de conclusie van deze masterproef.

Hoofdstuk 2 Gerelateerd werk Dit hoofdstuk geeft een kort overzicht van de nodige technische achtergrond en het gerelateerd werk voor deze masterproef. Eerst wordt in Paragraaf 2.1 de grafische processor besproken. Vervolgens behandelt Paragraaf 2.2 de dualiteit van prestatie en productiviteit. Daarna wordt de hoogniveau programmeertaal Julia voorgesteld in Paragraaf 2.3. Tot slot bespreekt Paragraaf 2.4 de tracetransformatie, welke zal gebruikt worden voor het evalueren van dit onderzoek.

2.1

Grafische processor

De GPU vindt zijn origine in grafische toepassingen: de ontwikkeling van grafische interfaces voor besturingssystemen en de popularisatie van 3D graphics voor allerhande toepassingen, creëerde een markt voor chips die gespecialiseerd zijn in het berekenen van beelden. Dankzij de grote populariteit van de eerste first person shoorters zoals Doom en Quake, ontstond een zoektocht naar het creëren van alsmaar realistischere omgevingen voor computerspellen. Zo evolueerde de GPU al snel tot zeer krachtige hardware. Het programmeren van deze hardware verliep echter door middel van een geparametriseerde pijplijn: de programmeur moest deze pijplijn enkel configureren (bv. de correcte lichtbron instellen) en van de juiste data voorzien (punten in een 3D ruimte). Later evolueerde deze pijplijn naar een programmeerbare pijplijn: voor bepaalde stappen in de pijplijn kon code geschreven worden die door de GPU uitgevoerd werd. Dit leidde ertoe dat wetenschappers een GPU begonnen te gebruiken voor algemene berekeningen los van enige grafische context. De enige manier om een GPU te programmeren was echter nog steeds via een grafische API, wat een aantal problemen met zich 4

2.1. GRAFISCHE PROCESSOR

5

meebracht. Ten eerste was kennis van deze uitgebreide API en de specifieke programmeertalen voor die stappen vereist en moesten algoritmes en probleemstellingen zodanig omschreven worden dat ze via deze API uit te voeren waren. Daarnaast was de GPU enkel bedoeld voor het uitvoeren van die grafische pijplijn en de hardware bevatte dus verschillende restricties voor het uitvoeren van algemene berekeningen. Hoewel initiële experimenten veelbelovende resultaten opleverden, introduceerde dit programmeermodel teveel restricties om algemeen toegepast te worden [5, 6, 7]. Als antwoord op deze beperkingen kwam NVIDIA met CUDA (Compute Unified Device Architecture), een parallel rekenplatform en programmeermodel [8]. Enerzijds bevat het de CUDA architectuur, een gewijzigde GPU architectuur die nieuwe componenten bevat. Het doel hiervan is het verwijderen van de tot dan toe aanwezige hardware restricties voor het uitvoeren van generieke algoritmes op de GPU, zonder de grafische capaciteiten in het gedrang te brengen. Anderzijds werd ook de taal CUDA C ontwikkeld, een extensie van standaard C, die het mogelijk maakt om los van de grafische API of pijplijn algoritmes voor de GPU te programmeren. Later werd deze extensie ook doorgevoerd op andere talen, zoals C++ en Fortran. CUDA is echter niet het enige platform dat het mogelijk maakt om een GPU voor algemene berekeningen te gebruiken. Een bekend alternatief, ontwikkeld door Apple en gestandaardiseerd door Khronos Group, is OpenCL (Open Computing Language). Dit is een raamwerk dat het mogelijk maakt om zogeheten kernels, geschreven in een taal gebaseerd op C99, uit te voeren op verschillende types hardware, bijvoorbeeld een CPU of GPU. Een beperking van dit raamwerk met betrekking tot het onderzoek in deze masterproef is het feit dat de OpenCL runtime de broncode van een kernel verwacht en geen gecompileerde code. Aangezien de bedoeling is om Julia code uit te voeren op de GPU, impliceert deze vereiste dat een vertaler geschreven wordt die Julia code naar OpenCL code vertaalt. Een dergelijke stap is bij gebruik van CUDA niet nodig, aangezien de CUDA runtime wel gecompileerde code aanvaardt.

2.1.1

Het programmeren van een GPU met CUDA

GPGPU programmeren aan de hand van CUDA gebeurt op basis van een heterogene architectuur. De opstelling bestaat uit een device, de GPU met bijhorend geheugen, dat de parallelle code uitvoert en een host, de CPU met bijhorend geheugen, dat seriële code uitvoert en het device beheert. De interactie tussen beide gebeurt ruwweg in 3 stappen: De invoerwaarden worden van het host geheugen naar het device geheugen gekopieerd De parallelle code wordt in de device geladen en uitgevoerd

6

HOOFDSTUK 2. GERELATEERD WERK De resultaatwaarden worden terug naar het host geheugen gekopieerd

Beschrijven welke berekeningen de GPU dient uit te voeren gebeurt aan de hand van een kernel. Dit is een C functie die omschrijft wat een draad op de GPU moet uitvoeren en wordt gekenmerkt door de global declaratie. Op de GPU worden vervolgens een aantal draden opgestart waarbij elke draad deze kernel uitvoert. De <<<...>>> syntaxis wordt gebruikt om op te geven hoeveel draden voor een bepaalde kernel opgestart moeten worden. Op software niveau worden de draden georganiseerd in een hiërarchische structuur: draden worden gegroepeerd in een blok en blokken worden gegroepeerd in een rooster. Dit is te zien in Figuur 2.1. Voor elke kernel wordt een één-, twee- of driedimensionaal rooster geconfigureerd. Het aantal blokken per roosterdimensie wordt ingesteld aan de hand van een CUDA C datatype: dim3 blokkenPerRooster( mx, my, mz );. Voor een één- of tweedimensionale structuur worden dan slechts één of twee argumenten respectievelijk opgegeven. Blokken bevatten eveneens een één-, twee- of driedimensionale structuur waarbij hetzelfde datatype gebruikt wordt om het aantal draden per blokdimensie te configureren: dim3 dradenPerBlok( nx, ny, nz );. Deze driedimensionale, hiërarchische organisatie maakt het mogelijk om op eenvoudige wijze één-, twee- of driedimensionale datastructuren te indexeren. Om in een kernel te kunnen bepalen welke draad binnen welk blok deze specifieke instantie uitvoert, zijn er in CUDA C een aantal variabelen gedefinieerd. De index van een draad kan opgevraagd worden aan de hand van de threadIdx variabele. Aangezien draden in een blok driedimensionaal georganiseerd worden, is deze variabele een vectortype met 3 elementen: threadIdx.[x|y|z]. Zo kunnen draden ook op driedimensionale wijze ge¨ındexeerd worden. Indexering van blokken binnen een rooster gebeurt eveneens op driedimensionale wijze aan de hand van een ingebouwd vectortype: blockIdx.[x|y|z]. Het aantal draden in elke dimensie van een blok en het aantal blokken in elke dimensie van een rooster kan opgevraagd worden door middel van de ingebouwde variabelen blockDim.[x|y|z] en gridDim.[x|y|z]. Een voorbeeld van een kernel en het instellen van de blok- en rooster configuratie is te zien in Codevoorbeeld 2.1. Hier worden twee vectoren van elk 72 elementen opgeteld waarbij elke draad op de GPU één paar elementen optelt. Er worden 2 × 3 blokken in het rooster en 3 × 4 draden per blok gealloceerd. Op basis van de draaden blokindex kan op eenvoudige wijze een uniek rangnummer voor elke draad berekend worden, zoals te zien in de eerste regel van de kernelfunctie in Codevoorbeeld 2.1.


Figuur 2.1: Voorbeeld van een 3 × 4 blokconfiguratie en een 2 × 3 roosterconfiguratie

__global__ void vadd(float* a, float* b, float* c) { int i = threadIdx.x + blockIdx.x * blockDim.x; c[i] = a[i] + b[i]; } int main() { ... dim3 blokkenPerRooster( 2, 3 ); dim3 dradenPerBlok( 3, 4 ); vadd<<< blokkenPerRooster, dradenPerBlok >>>( a, b, c ); ... } Codevoorbeeld 2.1: Voorbeeld van een GPU kernel, geschreven in CUDA C

7

8

HOOFDSTUK 2. GERELATEERD WERK

2.1.2

Compilatieproces en interactie met de GPU

Zoals vermeld in Paragraaf 2.1.1, bestaat de code van een programma dat door middel van CUDA gebruik maakt van de GPU, uit twee delen: host en device code. Compilatie van een dergelijk programma gebeurt aan de hand van de NVIDIA CUDA compiler. Host en device code worden door deze compiler apart verwerkt: Voor het compileren van host code zal de CUDA compiler enkel de <<<...>>> syntaxis vervangen door de CUDA runtime API functies die nodig zijn voor het inladen en opstarten van een kernel. De resulterende code kan dan verder gecompileerd worden aan de hand van een standaard C compiler, zoals GCC. Device code wordt door de CUDA compiler gecompileerd naar PTX instructies. De PTX (Parallel Thread Execution) ISA is een virtuele ISA en wordt door de CUDA driver tijdens uitvoering vertaald naar machine-instructies specifiek voor de GPU waarop de code wordt uitgevoerd. Het belangrijkste doel van deze ISA is GPU code porteerbaar te maken over verschillende GPU generaties.

De host is verantwoordelijk voor het beheer van het device. Hiervoor wordt in de host code gebruik gemaakt van de CUDA runtime API. De CUDA runtime, die de runtime API implementeert, is gebouwd op basis van de CUDA driver API. Dit is te zien in Figuur 2.2. Deze API biedt meer controle over laag niveau concepten zoals contexten, het equivalent van een CPU proces, en modules, het equivalent van dynamisch ingeladen bibliotheken. In dit geval moeten kernels nog steeds gecompileerd worden aan de hand van de CUDA compiler, maar de host code niet meer. Elke taal of omgeving kan gebruikt worden voor de interactie met de GPU, zolang deze de CUDA driver API kan aanspreken. De meeste toepassingen maken niet rechtstreeks gebruik van deze API, omdat dergelijk niveau van controle niet nodig is. De CUDA runtime API biedt dan meer gemak door contexten modulebeheer impliciet uit te voeren. In deze masterproef zal echter wel gebruik gemaakt worden van de CUDA driver API, aangezien de interactie met de GPU vanuit Julia zal verlopen.

2.1.3

GPU architectuur

Een schematische weergave van de belangrijkste componenten in de CUDA architectuur staat weergegeven in Figuur 2.3. Hier is te zien dat een GPU is opgebouwd uit meerdere multiprocessors (SM’s), waarvan het specifiek aantal afhankelijk is van generatie en model van de GPU. Elke multiprocessor bevat een aantal CUDA cores. Dit aantal is steeds een veelvoud van 32 en de controle over deze eenheden verloopt in groepen van 32. Met andere woorden, alle cores binnen een groep van 32 voeren dezelfde instructie uit,


9

Figuur 2.2: Organisatie van de CUDA runtime en de CUDA driver

mogelijks op verschillende data. Elke CUDA core bevat een volledig gepijplijnde ALU (Arithmetic Logic Unit) en FPU (Floating Point Unit). Naast deze cores zijn ook een aantal SFU’s (Special Function Unit) aanwezig. Deze behandelen meer complexe berekeningen, zoals sinus of cosinus. Ook SFU’s worden per 32 gecontroleerd. Tenslotte bevat de multiprocessor ook Load/Store eenheden voor het berekenen van bron en doeladressen. Wanneer de CPU een kernel opstart op de GPU, worden de blokken in het rooster van deze kernel verdeeld over de verschillende multiprocessors. Dit gebeurt op een schaalbare wijze, zoals te zien in Figuur 2.4. Hoeveel blokken tegelijk op een multiprocessor uitgevoerd kunnen worden is afhankelijk van de hardware die de draden in dat blok nodig hebben, zoals het aantal registers. Wanneer op een multiprocessor niet voldoende hardware aanwezig is om een blok van een kernel uit te voeren, faalt het opstarten van de kernel. Indien mogelijk worden meerdere blokken tegelijk op eenzelfde multiprocessor uitgevoerd. Wanneer een blok afgewerkt is, wordt een nieuw blok aan de multiprocessor toegevoegd. Dit levert een verzameling draden op die tegelijk op de multiprocessor aanwezig zijn. Het beheer van deze draden gebeurt in groepen van 32, genaamd warps. De draden binnen een warp worden uitgevoerd uit volgens een SIMT-principe (Single Instruction, Multiple Thread): de draden voeren dezelfde instructie uit, maar op verschillende data. Elke klokcyclus wordt door elke warp selector op de multiprocessor een uitvoerbare warp geselecteerd, waarvan dan, afhankelijk van de generatie van de architectuur, tot 2 onafhankelijke

10


Figuur 2.3: Schematische weergave van de verschillende componenten op een GPU

instructies opgehaald en uitgevoerd worden. De uitvoeringscontext van een warp bevat o.a. de programmatellers en registers van die warp en wordt op de multiprocessor bijgehouden gedurende de levensduur van de warp. Het wisselen tussen uitvoeringscontexten brengt dus geen kost met zich mee, wat wil zeggen dat elke klokcyclus verschillende warps geactiveerd kunnen worden. Dit laat toe om vertraging binnen een warp zoveel mogelijk te verbergen. Een multiprocessor bevat één registerbestand dat verdeeld wordt onder de verschillende warps en gedeeld geheugen dat verdeeld wordt onder de verschillende blokken op die multiprocessor. Het aantal blokken en warps van een kernel dat tegelijk door een multiprocessor verwerkt kan worden, is dus afhankelijk van het aantal registers en de hoeveelheid gedeeld geheugen dat voor die kernel nodig is. Twee draden die tot dezelfde warp behoren zullen steeds fysiek in parallel worden uitgevoerd. De relatieve ordening van instructies is dus steeds gewaarborgd: bij het uitvoeren van een instructie zijn gegarandeerd alle voorgaande instructies door alle draden van die warp uitgevoerd. Wanneer de twee draden echter tot verschillende warps behoren, geldt deze garantie niet langer. Indien er afhankelijkheden zijn tussen draden, kan dit problemen geven. Een voorbeeld hiervan is het uitlezen van een waarde uit het geheugen die door een andere draad werd weggeschreven. Op het moment dat de draad de waarde probeert uit te lezen, is deze niet noodzakelijk aanwezig in het geheugen. Dergelijke afhankelijkheden tussen draden kunnen gewaarborgd worden door gebruik te maken van de ingebouwde syncthreads() functie. Deze functie geldt als een barrière in de code waar alle draden moeten aankomen vooraleer de uitvoering van de draden hervat mag worden. De GPU bevat verschillende types geheugen die voor een draad toegankelijk zijn om data naar te schrijven en van te lezen. Elke draad beschikt over lokaal geheugen dat

2.2. PRODUCTIVITEIT VERSUS PRESTATIE

11

Figuur 2.4: Blokken van een kernel worden op dynamische wijze verdeeld over de aanwezig multiprocessors

enkel toegankelijk is voor die draad. Daarnaast heeft elk blok een hoeveelheid gedeeld geheugen dat toegankelijk is voor elke draad in dat blok. Tenslotte hebben alle draden ook nog toegang tot het globaal geheugen. Dit globaal geheugen is gescheiden van het geheugen van de CPU. Aangezien de GPU geen toegang heeft tot het geheugen van de CPU en dus ook geen data van dit geheugen naar het globaal GPU geheugen kan sturen, valt het beheer van het globaal geheugen van de GPU onder de verantwoordelijkheid van de CPU. Dit houdt het alloceren en vrijgeven van geheugen in, alsook het transport van data tussen CPU en GPU geheugen. Een zeer recentelijke toevoeging aan CUDA is het principe van unified memory. Dit biedt aan de gebruiker de illusie dat de CPU en GPU gebruik maken van hetzelfde geheugen. De gebruiker hoeft in de host code enkel nog geheugen te alloceren en vrij te geven. Intern is het geheugen echter nog steeds gescheiden en is de CUDA runtime verantwoordelijk voor het transport van data tussen beide.

2.2

Productiviteit versus prestatie

Ondanks de hoge prestaties die gehaald kunnen worden met CUDA op een GPU, is het gebruik ervan complex, tijdrovend of foutgevoelig. Daarom proberen verschillende

12


onderzoeksprojecten het programmeren van een GPU productiever te maken, bijvoorbeeld [9]. Meer in het algemeen zijn laagniveau talen niet productief omdat de programmeur geconfronteerd wordt met te veel technische details. In wetenschappelijke en wiskundige context wordt daarom steeds meer gekozen om te werken in een hoogniveau, ge¨ınterpreteerde programmeertaal. Bekende voorbeelden hiervan zijn MATLAB, R [10] en SciPy [11]. Het succes van MATLAB en andere wetenschappelijke hoogniveau talen, is te danken aan het gemak en productiviteit dat ze bieden ten opzichte van laagniveau talen zoals C of Fortran. Voorbeelden hiervan zijn een eenvoudigere syntaxis, integratie met een omgeving voor simulatie en visualisatie, de mogelijkheid tot het interactief uitvoeren van commando’s, een uitgebreide verzameling numerieke functies die zeer eenvoudig gebruikt kunnen worden, etc. Een nadeel van deze hoogniveau talen is dat ze meestal niet in staat zijn om dezelfde prestaties te behalen als bij een laagniveau taal. Om dit verlies aan prestatie tegen te gaan, werd in veel hoogniveau talen het principe van een tweestaps ontwikkeling van algoritmes toegepast: de logica van het algoritme wordt in een hoogniveau taal geschreven, terwijl voor het echte rekenwerk een laagniveau taal zoals C of Fortran wordt gebruikt. Een voorbeeld hiervan is NumPy [12], een Python bibliotheek die een datastructuur aanbiedt om n-dimensionale rijen efficiënt voor te stellen en te bewerken. Deze bibliotheek wordt gebruikt in Python, terwijl ze grotendeels bestaat uit geoptimaliseerde C code. Hoewel tweestaps ontwikkeling in bepaalde gevallen een goede strategie is, heeft het toch een aantal nadelen. Zo kan er bij het implementeren van de rekenintensieve delen niet geprofiteerd worden van de productiviteit die een hoogniveau taal biedt. Ook moet er kennis zijn van beide talen en introduceert de interactie ertussen vaak een significante meerkost. Analyses en optimalisaties kunnen ook niet meer op het algoritme als geheel uitgevoerd worden. Een alternatieve oplossing voor de tweestaps ontwikkelingsstrategie is de prestatie van een bestaande taal zelf te verbeteren. Een bekend voorbeeld hiervan dat reeds goede resultaten heeft bereikt, is het Python compiler raamwerk PyPy [13]. Het grote voordeel van deze methode is dat vrijwel alle bestaande code in die taal van deze versnellingen kan profiteren. In de praktijk echter hebben dergelijke projecten de nood aan een tweestaps ontwikkelingsstrategie nog niet geëlimineerd. In sommige talen heeft men, ermee rekening houdend dat de taal op basis van een interpreter zou ge¨ımplementeerd worden, zelfs ontwerpbeslissingen genomen die het genereren van efficiënte code hinderen [14]. Julia is een hoogniveau programmeertaal die de tekortkomingen van beide voorgaande methodes aanpakt. De taal en implementatie ervan zijn van in het begin opgebouwd om zoveel mogelijk gebruik te maken van technieken om dynamische talen efficiënt uit te voeren. Ze biedt de productiviteit van een hoogniveau programmeertaal, terwijl de prestaties vergelijkbaar zijn met die van statisch gecompileerde talen. Paragraaf 2.3 komt hier uitvoerig op terug.

2.3. JULIA

13

Een typisch kenmerk van hogere programmeertalen, zoals Julia of MATLAB, is dat ze de productiviteit verhogen door de nodige abstracties aan te bieden aan de programmeur. Het gebruik van abstracties wordt in verschillende onderzoeksprojecten behandeld [15, 16]. Vooral in de context van parallelle berekeningen (multi-core, GPU, etc.) wordt dit beschouwd als dé strategie om de productiviteit van softwareontwikkeling te verbeteren. Deze masterproef bouwt hierop verder door abstracties aan te bieden in Julia om het gebruik van een GPU productiever te maken. In [17] probeert men in de programmeertaal Rust hetzelfde doel te bereiken op een gelijkaardige manier. Deze masterproef verschilt echter op twee belangrijke punten. Ten eerste wordt Julia dynamisch gecompileerd waardoor runtime informatie kan gebruikt worden om de abstractie te optimaliseren; daarentegen wordt Rust statisch gecompileerd en is dergelijke informatie niet ter beschikking. Ten tweede is de programmeur in de Rust-oplossing verantwoordelijk voor de administratieve taken zoals geheugenbeheer en het laden van de code in de GPU, terwijl in deze masterproef voorzieningen genomen zijn om dergelijke taken transparant te maken.

2.3 2.3.1

Julia De taal

Julia [4] is een hoogniveau, dynamische programmeertaal die ontwikkeld werd aan MIT en vooral gericht is op het gebruik in wetenschappelijke en wiskundige context. De taal is ontwikkeld met als doel productiviteit én prestatie. Het succes van dit project is te zien in Figuur 2.5. De belangrijkste aspecten die leiden tot deze goede prestatie zijn: Julia code wordt niet ge¨ınterpreteerd, maar tijdens uitvoering gecompileerd. Hiervoor wordt gebruik gemaakt van de LLVM JIT compiler. LLVM [18] (Low Level Virtual Machine) is een modulaire compilatie-infrastructuur die een taal- en platformonafhankelijke intermediaire representatie van code aanbiedt waarop analyses en optimalisaties uitgevoerd kunnen worden. Code wordt op het moment van uitvoeren geoptimaliseerd op basis van typeinformatie die op twee manieren verzameld wordt. Enerzijds is het principe van multiple dispatch een kernonderdeel van Julia. Dit maakt het voor de gebruiker mogelijk om, indien gewenst, het gedrag van een functie aan te passen op basis van typedeclaraties van de argumenten. Op deze manier geeft de programmeur de compiler al type-informatie nog voor de functie uitgevoerd wordt. Anderzijds is het type van waarden tijdens uitvoering gekend en kan ook deze informatie gebruikt worden voor optimalisatie.

14


103 102 101 100

fib

parse int Fortran

Julia

quicksort Python

mandel

MATLAB

Figuur 2.5: Benchmark uitvoeringstijden relatief ten opzichte van C

Er bestaan reeds veel hoogefficiënte bibliotheken, zoals BLAS en LAPACK. Om niet steeds het wiel opnieuw uit te vinden, kunnen deze bibliotheken met quasi geen meerkost vanuit Julia aangesproken worden.

2.3.2

Werking van de Julia compiler

De compilatie van Julia code verloopt in verschillende stappen. Eerst wordt de Julia code verwerkt door een parser die twee taken uitvoert. Enerzijds wordt de syntaxis van de Julia code omgezet naar een meer eenvoudige representatie die het gemakkelijker maakt om analyses uit te voeren en code te genereren. Voorbeelden hiervan zijn het uit elkaar halen van bepaalde geneste expressies en het verwijderen van “syntactische suiker”. Tijdens deze omzetting worden ook bepaalde variabele en closure analyses uitgevoerd. Anderzijds genereert de parser een abstracte syntaxisboom (AST). Aangezien de parser in FemtoLisp ge¨ımplementeerd is, is deze AST een s-expressie, een populaire notatie voor boomgestructureerde data in Lisp. De overige stappen van de Julia compilatie zijn in C ge¨ımplementeerd. Bijgevolg wordt de FemtoLisp AST omgezet naar een in C ge¨ımplementeerde AST. In geval van een toplevel expressie wordt deze AST verder gecompileerd en uitgevoerd. Beschrijft de AST echter een functie, dan wordt de AST omgezet in een methodedefinitie en toegevoegd aan een methodetabel. Voor elke functie wordt een methodetabel aangemaakt die de verschillende versies van die functie bevat. De verschillende versies worden gekenmerkt door een verschil in aantal en/of type van de argumenten van de functie.

2.3. JULIA

15

Wanneer een functie dan opgeroepen wordt, dient de AST ervan verder gecompileerd te worden. Op basis van het type van de concrete argumentwaarden, wordt de juiste versie opgezocht in de methodetabel. Op basis van de bijhorende AST worden LLVM IR instructies gegenereerd. Hierbij wordt elke knoop in de AST omgezet in de nodige instructies. Dit gebeurt aan de hand van een IRBuilder die een uniforme API aanbiedt voor het aanmaken van instructies en het toevoegen van deze instructies in basisblokken. Op de gegenereerde LLVM IR kunnen vervolgens optimalisaties uitgevoerd worden. Dit gebeurt aan de hand van een PassManager en is te vergelijken met de --On optie van bijvoorbeeld GCC. De gegenereerde LLVM IR instructies worden tenslotte gecached in de methodetabel om te vermijden dat bij het opnieuw oproepen van een functie, wederom deze compilatiestap uitgevoerd zou moeten worden.

In de laatste stap moeten de hoogniveau LLVM IR instructies omgezet worden in machinecode voor een specifiek platform. Hiervoor wordt een gerichte acyclische graaf opgesteld die voor elke LLVM IR instructie een knoop bevat. De belangrijkste informatie in zo’n knoop is de opcode en de argumenten van de operatie van die knoop. Pijlen tussen knopen drukken data- en controleverloop afhankelijkheden uit. De opgestelde graaf biedt een abstracte representatie van code waarop instructieselectie kan uitgevoerd worden met automatische technieken. Aan deze graaf worden vervolgens zoveel mogelijk platformspecifieke details toegevoegd, om laagniveau optimalisaties uit te kunnen voeren. Op de resulterende graaf kan echter nog geen instructieselectie toegepast worden aangezien de graaf nog niet legaal is: niet alle gebruikte types of operaties worden door het platform ondersteund en moeten dus aangepast worden, wat in de legalisatiefase gebeurt. Na deze fase wordt de graaf opnieuw geoptimaliseerd om de door de legalisatiefase ge¨ıntroduceerde inefficiënties weg te werken. Vervolgens wordt op basis van de legale graaf instructieselectie uitgevoerd. Hierbij worden operaties of combinaties ervan uit de legale graaf omgezet in concrete machine-instructies die samen eveneens een gerichte acyclische graaf vormen. In de planningsfase wordt deze graaf omgezet in een lineaire sequentie van machine-instructies. Bepaalde heuristieken worden hier toegepast om snellere code te verkrijgen dan het geval zou zijn bij het simpelweg topologisch sorteren van de instructies uit de graaf. Tenslotte dient ook nog registerallocatie uitgevoerd te worden. De machine-instructies maken op dit punt immers nog gebruik van een oneindig aantal virtuele registers, behalve die instructies die een specifiek register vereisen. De resulterende sequentie van machine-instructies kan dan weggeschreven worden naar een objectbestand.

16

2.4


Tracetransformatie

Om het werk van deze masterproef te evalueren, wordt de oplossing gebruikt voor een concrete toepassing: de tracetransformatie [19]. Dit is een transformatie die eigenschappen uit figuren probeert te extraheren, waarbij deze eigenschappen invariant zijn onder bepaalde beeldtransformaties, zoals rotatie, translatie en schaling. De tracetransformatie wordt in verschillende toepassingen voor objectherkenning en -identificatie gebruikt, zoals gezichtsherkenning [20] en het schatten van parameters van geometrische transformaties [21]. Concreet worden de eigenschappen van een figuur in verschillende stappen berekend, wat visueel weergegeven is in Figuur 2.6. Eerst wordt een figuur over een hoek geroteerd. Op de kolommen van de resulterende figuur wordt dan een T functionaal toegepast, die voor elke kolom een scalaire waarde berekent. Voor elke rotatie wordt zo een vector van scalairen verkregen, één voor elke kolom. Al deze vectoren samen vormen een tweedimensionale matrix, het sinogram. Deze kan vervolgens verder verwerkt worden door een P functionaal toe te passen op de verschillende vectoren, om zo per vector opnieuw een scalair te verkrijgen. In een laatste stap kan deze vector geaggregeerd worden tot een enkele scalair door toepassing van een Φ functionaal. Hoewel de robuustheid van de tracetransformatie werd aangetoond [19], wordt ze nog maar weinig toegepast in de praktijk omwille van de hoge computationele complexiteit [22]. Aangezien verschillende delen in de tracetransformatie mogelijkheden tot parallelle uitvoer vertonen, kan het gebruik van een GPU een hoge prestatiewinst opleveren. In veel situaties waar de tracetransformatie wordt toegepast, ontbreekt echter de kennis of de tijd om de implementatie opnieuw uit te voeren. Deze toepassing leent zich dus zeer goed tot het evalueren van de oplossing van dit onderzoek.

2.4. TRACETRANSFORMATIE

17

Figuur 2.6: Schematische weergave van de verschillende stappen in de tracetransformatie

Hoofdstuk 3 Compilatie van Julia naar PTX Dankzij de populariteit van het gebruik van GPU’s in wetenschappelijke context, zijn voor veel programmeertalen pakketten geschreven die via de CUDA driver API rechtstreeks toegang bieden tot de GPU. Ook in Julia bestaat zo’n pakket: CUDA.jl [23]. Het gebruik van dit pakket verloopt in meerdere stappen. Ten eerste moet een GPU kernel geschreven worden in CUDA C. In Codevoorbeeld 3.1 is een kernel te zien die twee matrices optelt. Vervolgens moet de programmeur deze code compileren naar een PTX module door middel van de CUDA compiler. Daarna kan de geschreven kernel uitgevoerd worden vanuit Julia code, wat te zien is in Codevoorbeeld 3.2. Hoewel dit pakket veel nuttige functionaliteit aanbiedt, vereist het wel de nodige kennis: de programmeur is immers zelf verantwoordelijk voor het compileren van de kernels, het inladen van modules, alloceren en transporteren van geheugen naar de GPU en terug enz. Om deze negatieve impact op de productiviteit tegen te gaan, werd in deze masterproef een infrastructuur ontwikkeld die aan volgende vereisten voldoet: Kernels moeten in Julia geschreven kunnen worden en de compilatie ervan dient automatisch te gebeuren. Hierdoor hoeft de programmeur zelf geen CUDA C te kennen of te schrijven.

__global__ void vadd(const float *a, const float *b, float *c) { int i = threadIdx.x + blockIdx.x * blockDim.x; c[i] = a[i] + b[i]; } Codevoorbeeld 3.1: Voorbeeld van een GPU kernel, geschreven in CUDA C

18

19

using CUDA # PTX code inladen en functie opzoeken md = CuModule("vadd.ptx") vadd = CuFunction(md, "vadd") # genereer de input rijen a = round(rand(Float32, (3, 4)) * 100) b = round(rand(Float32, (3, 4)) * 100) # transporteer ze naar het GPU geheugen ga = CuArray(a) gb = CuArray(b) # alloceer op het GPU geheugen een rij # die de resultaatwaarden zal bevatten gc = CuArray(Float32, (3, 4)) launch(vadd, 3, 4, (ga, gb, gc)) # transporteer de resultaten terug naar het host geheugen c = to_host(gc) free(ga) free(gb) free(gc) unload(md) Codevoorbeeld 3.2: Gebruik van CUDA.jl in Julia

20

HOOFDSTUK 3. COMPILATIE VAN JULIA NAAR PTX Administratieve taken die gepaard gaan met het uitvoeren van een kernel op de GPU moeten geautomatiseerd worden. Dit omvat o.a. het inladen van de module op de GPU, het alloceren en vrijgeven van geheugen op de GPU en het transport van geheugen naar de GPU en terug.

Om de integratie met Julia zo vlot mogelijk te laten verlopen en de productiviteit van de programmeur niet te hinderen, was het belangrijk om zoveel mogelijk de bestaande taalabstracties aanwezig in Julia te behouden en ze zo goed mogelijk te mappen op GPU functionaliteit. In dit hoofdstuk wordt de infrastructuur die in deze masterproef ontwikkeld werd, voorgesteld. In Paragraaf 3.1 wordt de belangrijkste functionaliteit ervan vermeld. Paragraaf 3.2 beschrijft beperkingen van de huidige Julia compiler, rekening houdend met de functionaliteit die nodig is voor deze infrastructuur. Tenslotte wordt in Paragraaf 3.3 besproken hoe deze infrastructuur intern ge¨ımplementeerd werd en welke aanpassingen gemaakt zijn aan de Julia compiler.

3.1

API in Julia

Het gebruik van onze infrastructuur is vergelijkbaar met die van CUDA.jl. Eerst dient in Julia een kernel functie geschreven te worden die bestemd is voor de GPU. Deze kernel functie is zeer gelijkaardig aan een gewone Julia functie, op één verschil na. In Julia wordt er automatisch van uit gegaan dat de waarde van de laatste expressie ook het resultaat van de functie is, tenzij expliciet anders vermeld met een return statement. Een host kernel op de GPU zal de berekende resultaten echter rechtstreeks naar het geheugen schrijven zonder deze expliciet terug te geven; een host kernel moet dus steeds een void functie zijn. Bijgevolg dient dit ook expliciet in de Julia functie vermeld te worden door een nothing waarde terug te geven. Aangezien er geen resultaten teruggeven worden, moeten de variabelen die de resultaten dienen te bevatten ook als argument in de functiedeclaratie opgegeven worden. Een voorbeeld van een in Julia geschreven kernel is te zien in Codevoorbeeld 3.3. Deze kernel functie moet zich bevinden in een Julia module waarvan de naam begint met ptx ; de reden hiervoor wordt besproken in Paragraaf 3.3.1. Om de functienaam ook buiten de module beschikbaar te maken, dient de functie geëxporteerd te worden. Voor het implementeren van kernel functies werd in deze masterproef volgende functionaliteit aan Julia toegevoegd: Er werd ondersteuning toegevoegd voor dynamisch gedeeld geheugen:

3.1. API IN JULIA

21

module __ptx__Example export vadd function vadd(a,b,c) i = threadId_x() + blockId_x() * numThreads_x() c[i] = a[i] + b[i] return nothing end end Codevoorbeeld 3.3: Voorbeeld van een GPU kernel, geschreven in Julia

– shmem = cuSharedMem(): declareer een variabele als dynamisch gedeeld geheugen. – setCuSharedMem(shmem, index, value): schrijf een waarde weg naar het voorheen gedefinieerd dynamisch gedeeld geheugen op een bepaalde positie. – getCuSharedMem(shmem, index): lees een waarde uit het voorheen gedefinieerd dynamisch gedeeld geheugen op een bepaalde positie. Bij het schrijven van kernels in CUDA C zijn een aantal intrinsics voorzien die informatie verschaffen over de uitvoerende context. Dezelfde intrinsics werden ook voorzien in Julia code:

– threadId x/ y/ z() voor het opvragen van de index van een draad in elke dimensie – numThreads x/ y/ z() voor het opvragen van het aantal draden per blok in elke dimensie – blockId x/ y/ z() voor het opvragen van de index van een blok in elke dimensie – numBlocks x/ y/ z() voor het opvragen van het aantal blokken in elke dimensie – sync threads() voor het synchroniseren van alle draden in eenzelfde blok Ook werd toegang tot een aantal wiskundige functies voorzien, namelijk sin(), cos() en floor(). De reden waarom de standaard wiskundige functies in Julia niet gebruikt kunnen worden, wordt uitgelegd in Paragraaf 3.2.6.

22

HOOFDSTUK 3. COMPILATIE VAN JULIA NAAR PTX

using __ptx__Example a = round(rand(Float32, (3, 4)) * 100) b = round(rand(Float32, (3, 4)) * 100) c = Array(Float32, (3,4)) @cuda (__ptx__Example, 3, 4) vadd(CuIn(a),CuIn(b),CuOut(c)) Codevoorbeeld 3.4: Voorbeeld van het gebruik van de @cuda macro

Om de in Julia geschreven kernel uit te voeren werd een macro geschreven die ook de administratieve taken uitvoert. Dit is te zien in Codevoorbeeld 3.4. Deze macro vereist twee argumenten: Een set configuratieparameters: ten eerste de Julia module waartoe de kernel behoort, daarnaast ook nog de draaden blokconfiguratie en tenslotte ook het aantal bytes dynamisch gedeeld geheugen. De kerneloproep zelf met concreet ingevulde argumenten: zoals hiervoor vermeld, moeten zowel de input- als de outputargumenten opgegeven worden. Om aan de macro duidelijk te maken welke argumenten als input en/of output gebruikt zullen worden, dienen de argumenten in een wrapper geplaatst te worden. De reden hiervoor wordt verklaard in Paragraaf 3.3.5.

3.2

Beperkingen van de Julia compiler

Julia is gebouwd op basis van de LLVM compilerinfrastructuur en de Julia frontend genereert dus LLVM IR die door bijvoorbeeld de x86 backend omgezet wordt in machinecode. Hoewel er in LLVM ook een backend voorzien is voor het compileren naar PTX instructies, nl. NVPTX, levert het simpelweg verwisselen van de backends geen geldige PTX code. Dit komt omdat de NVPTX backend niet de volledige LLVM IR ondersteunt. Zo ontbreekt bijvoorbeeld ondersteuning voor intrinsics die het optreden van een overflow bij bepaalde berekeningen herkennen. Anderzijds is de LLVM IR ook niet volledig target onafhankelijk. De semantiek van een fence instructie is bijvoorbeeld niet mogelijk op een GPU. Synchronisatie van draden dient daardoor te gebeuren aan de hand van de daarvoor voorziene PTX intrinsics. De voor de GPU bestemde code moet dus apart beheerd en verwerkt worden. Bijgevolg kon de Julia frontend wel gebruikt worden voor het genereren van NVPTX compatibele IR, op voorwaarde dat in de codegeneratiestap een aantal wijzigingen werden aangebracht voor volgende beperkingen.

3.2. BEPERKINGEN VAN DE JULIA COMPILER

3.2.1

23

Gescheiden codebeheer

In LLVM wordt gebruik gemaakt van modules om eenheden van code aan te duiden, waarbij die code als een geheel verwerkt dient te worden. Een LLVM module bevat zaken zoals globale variabelen, functiedefinities en -implementaties. Aangezien de voor de GPU bestemde code anders verwerkt moet worden dan de gewone Julia code, dient deze ondergebracht te worden in een aparte LLVM module. Deze module wordt door de NVPTX backend gecompileerd naar PTX code, welke dan aan de hand van de CUDA driver API in de GPU geladen kan worden. Belangrijk is dus dat de gebruiker op een intu¨ıtieve manier moet kunnen aangeven welke functies als GPU kernels ge¨ınterpreteerd moeten worden.

3.2.2

Onderscheid tussen host en device functies

In PTX wordt een onderscheid gemaakt tussen host en device functies. Hierbij vormt een host functie het toegangspunt van een kernel en kan door de CPU aan de hand van de CUDA driver opgestart worden. Device functies zijn niet rechtstreeks toegankelijk voor de CPU en kunnen enkel opgeroepen worden vanuit GPU code. Om dit onderscheid te kunnen maken vereist de NVPTX backend dat de LLVM IR functies die als host functie bedoeld zijn, van een extra annotatie worden voorzien.

3.2.3

NVPTX Intrinsics

Zoals beschreven in Paragraaf 2.1, biedt CUDA C een aantal ingebouwde variabelen en functies aan. Deze worden in de NVPTX backend aangeboden als intrinsics, welke in LLVM IR te gebruiken zijn. Aangezien de programmeur deze variabelen en functies moet kunnen gebruiken bij het schrijven van een kernel, moeten deze intrinsics aangeboden worden in Julia.

3.2.4

Alternatieve telling

Net als bij sommige andere hoogniveau talen, zoals MATLAB, start de indexering van rijen en tellingen in het algemeen bij Julia vanaf 1, waar dit in CUDA C begint bij 0. Omdat de programmeur in Julia natuurlijk gewend is aan eerstgenoemde telling, dient hier ook rekening mee gehouden te worden. Dit geldt niet alleen bij het indexeren van rijen, maar ook bij het indexeren van draden en blokken. Wanneer de programmeur immers in een kernel de index van een draad of blok opvraagt, verwacht hij dat de eerste draad of het eerste blok index 1 heeft, terwijl de NVPTX intrinsics hiervoor 0 teruggeven.

24

3.2.5


Geheugenbeheer

Op een GPU zijn, zoals vermeld in Paragraaf 2.1.3, verschillende types geheugen aanwezig, die in LLVM IR gerepresenteerd worden door adresruimtes. Elke geheugenreferentie in LLVM IR dient dus voorzien te zijn van de addrspace annotatie, wat aanduidt tot welke adresruimte deze behoort. Dit gebeurt bijvoorbeeld bij een kernel die een rij als argument heeft. Die rij wordt door de host gealloceerd in het globaal geheugen van de GPU en dit moet expliciet in de LLVM code weergegeven worden. Naast globaal geheugen is er, zoals beschreven in Paragraaf 2.1.3, ook gedeeld geheugen. Dit geheugen is enkel toegankelijk voor alle draden binnen eenzelfde blok en kan gebruikt worden om draden binnen een blok met elkaar te laten communiceren. In geval van dynamisch gedeeld geheugen specificeert de programmeur bij het opstarten van de kernel hoeveel bytes hiervoor gealloceerd moeten worden. Zoals aangegeven in Paragraaf 3.1, moet de programmeur in Julia code dynamisch geheugen kunnen declareren en indexeren. Een tweede aspect van geheugenbeheer op een GPU is het alloceren van geheugen en de transfer van data tussen het CPU geheugen en het GPU geheugen. Zoals aangegeven in Paragraaf 3.1, moet de programmeur bij het oproepen van een kernel in Julia kunnen aangeven welke argumenten als input en/of ouput gebruikt worden door de kernel. Dit laat toe om onnodige geheugentransfers te vermijden. Tenslotte kan ook nog opgemerkt worden dat de programmeur gewend is dat Julia controleert wanneer een indexering van een rij buiten de grenzen van die rij valt. Dit gebeurt volledig automatisch, zonder dat de programmeur expliciet extra informatie dient te verschaffen. De datastructuur die Julia gebruikt voor een rij bevat namelijk, naast de rij zelf, extra informatie zoals het aantal dimensies en de grootte per dimensie. Op basis van deze informatie worden extra instructies gegenereerd die de nodige controles uitvoeren. Op de GPU is een rij echter simpelweg een geheugenreferentie naar de rij zelf, zonder de nodige extra informatie om dergelijke controles uit te voeren.

3.2.6

Linken met libdevice

Bij het schrijven van CPU code die gebruik maakt van wiskundige primitieven, zoals sinus en cosinus, dient gelinked te worden met libm, de standaardbibliotheek voor wiskundige functies. De code in deze bibliotheek is echter niet uitvoerbaar op de GPU. NVIDIA heeft daarom specifiek voor GPU’s de libdevice bibliotheek voorzien: een verzameling wiskundige en andere functies, geschreven in LLVM bitcode. De functies in deze bibliotheek zijn geoptimaliseerd voor specifieke GPU architecturen, dus na compilatie dient gelinked te worden met de juiste versie van libdevice, afhankelijk van de hardware waarop de code uitgevoerd zal worden. Deze stap wordt normaal gezien uitgevoerd door

3.3. AANPASSINGEN AAN DE JULIA COMPILER

25

de CUDA compiler. Voor het compileren van Julia kernels naar PTX instructies kan van deze compiler geen gebruik gemaakt worden en dient de Julia codegeneratie deze linker stap uit te voeren.

3.2.7

Een in Julia geschreven kernel opstarten

Wanneer een in Julia geschreven functie opgeroepen wordt aan de hand van de syntaxis (<argumenten...>), zal de compiler deze functie compileren en uitvoeren. De API uit Paragraaf 3.1 toont dat kernels als gewone functies geschreven worden. Deze kernels kunnen echter niet opgestart worden aan de hand van voornoemde syntaxis. De Julia compiler onderneemt namelijk geen acties om de code in de GPU te laden en andere administratieve taken uit te voeren. Bijgevolg moest extra functionaliteit voorzien worden om kernels vanuit Julia te kunnen opstarten.

3.3

Aanpassingen aan de Julia compiler

Paragraaf 3.2 gaf een overzicht van de beperkingen van de huidige Julia compiler voor wat betreft het gebruik van een GPU zoals gespecificeerd in Paragraaf 3.1. Hieronder wordt uitgelegd hoe deze beperkingen werden aangepakt.

3.3.1

Gescheiden codebeheer

De Julia compiler werd aangepast zodat deze nu, naast de standaard LLVM module, een tweede LLVM module bijhoudt die alle GPU kernels zal bevatten en door de NVPTX backend verwerkt zal worden. De programmeur moet dan op een intu¨ıtieve manier kunnen aangeven welke code is deze tweede LLVM module terecht dient te komen. Hiervoor wordt gebruik gemaakt van een Julia module. Dit is een constructie op code niveau die het mogelijk maakt om zonder naamconflicten variabelen, types en functies te definiëren en valt te vergelijken met het concept van namespaces in C++. De gebruiker kan dan alle GPU kernels in deze module schrijven. Om aan de Julia compiler aan te duiden welke Julia module de GPU kernels bevat, zou het interessant zijn mocht de gebruiker dit via een of ander attribuut kunnen specificeren, maar dit wordt voorlopig niet ondersteund door Julia. Bijgevolg is in dit eindwerk gekozen om de Julia module met GPU kernels te herkennen op basis van de naam. In het codegeneratieproces maakt Julia een datastructuur aan voor elke Julia module die gecompileerd wordt. Deze datastructuur werd uitgebreid zodat nu bij het aanmaken ervan wordt gekeken of de naam van de module begint met ptx . Of de module dan kernels bevat of niet wordt eveneens in die datastructuur bijgehouden. Bij de compilatie van een functie bevat de context van die functie

26


een referentie naar de datastructuur van de Julia module waartoe de functie behoort. Op die manier kan tijdens het codegeneratieproces van een functie eenvoudig gecontroleerd worden of de omvattende module kernels bevat en deze functie dus als kernel bedoeld is. De gegenereerde IR kan dan aangepast worden en in de voor de GPU bedoelde LLVM module gestoken worden.

3.3.2

Onderscheid tussen host en device functies

Om in de codegeneratie het onderscheid te kunnen maken tussen host en device functies, zou men hetzelfde principe kunnen hanteren als in CUDA C: de gebruiker duidt dit expliciet in de code aan. In Julia geeft dit echter twee problemen. Ten eerste zou zoiets de mogelijkheid vereisen om functieattributen te specificeren, wat momenteel niet ondersteund wordt door Julia. Ten tweede is het onderscheid tussen host en device functies een technisch detail dat enkel geldt voor GPU’s. Het doel van een hoogniveau taal zoals Julia, is juist dergelijke details verbergen, wat bij het het expliciet aanduiden van dit onderscheid niet gebeurt. Bijgevolg is gekozen om dit onderscheid op impliciete wijze te maken. Hiervoor wordt gebruik gemaakt van het feit dat functies die zichtbaar moeten zijn buiten de module waarin ze gedefinieerd zijn, geëxporteerd moeten worden. Functies die enkel binnen hun module gebruikt worden, dienen niet geëxporteerd te worden. Zoals beschreven in Paragraaf 3.3.1, worden GPU functies in een aparte Julia module gegroepeerd. Het opstarten van een GPU functie gebeurt in CPU code, buiten die module. Dit vereist dus dat die functie buiten zijn module zichtbaar is en dus geëxporteerd moet worden. Aangezien enkel host functies in CPU code opgestart worden, komt het exporteren van functies overeen met het onderscheid tussen host en device functies. Het exporteren van functies is echter een algemeen aspect van Julia, dus deze methode maakt dit onderscheid dus inderdaad op impliciete wijze. Deze oplossing heeft echter een beperking. Immers, indien er op termijn ook Julia bibliotheken geschreven worden die device functies aan programmeurs aanbieden, moeten de namen van deze functies buiten de module zichtbaar zijn, zodat de programmeur ze in zijn eigen GPU code gebruik van kan maken. Dit vereist dan dat deze functies geëxporteerd worden, terwijl het toch geen host functies zijn.

3.3.3

NVPTX Intrinsics

Een mogelijkheid om de NVPTX intrinsics in Julia aan te bieden is om hiervoor ook intrinsics in de Julia taal zelf te maken. Hoewel het toevoegen van intrinsics aan Julia mogelijk is, is dit expliciet moeilijk gemaakt omdat het niet de bedoeling is dat zomaar


27

nieuwe intrinsics aan de taal worden toegevoegd. De gesuggereerde methode is gebruik te maken van de voorlopig niet-officiële functionaliteit om rechtstreeks LLVM instructies te inlinen. Dit laat toe om de intrinsics rechtstreeks aan te spreken, zonder enige aanpassing aan de codegeneratie.

3.3.4

Alternatieve telling

Aangezien de telling in LLVM ook bij 0 start, was in Julia reeds de functionaliteit aanwezig die de juiste omzetting uitvoert. Dit houdt in dat bij elke rijtoegang een extra instructie gegenereerd wordt die 1 aftrekt van de index van deze toegang. Voor de rijtoegang in de GPU kernel kon dus deze functionaliteit herbruikt worden. Om ook de draaden blokindexering bij 1 te laten starten, wordt simpelweg 1 opgeteld bij de waarde die de intrinsics teruggeven.

3.3.5

Geheugenbeheer

Intern worden rijen in Julia niet voorgesteld door een geheugenreferentie, maar aan de hand van een aangepast datatype dat ook extra informatie over die rij bevat, zoals de lengte van de rij. Dit datatype wordt niet ondersteund door de GPU, een kernel verwacht dat de argumenten die een rij zijn doorgegeven worden als een eenvoudige geheugenreferentie naar het globaal geheugen. In de Julia codegeneratie wordt voor een kernel een gewijzigde LLVM functiedeclaratie gegenereerd. De argumenten van de functie worden overlopen en voor de argumenten die een rij zijn wordt als datatype een geheugenreferentie opgegeven in plaats van het specifiek Julia rijtype. Daarnaast wordt ook de addrspace(1) annotatie toegevoegd om aan te duiden dat de geheugenreferentie naar de globale adresruimte wijst. Zoals in Paragraaf 3.1 te zien is, wordt de ondersteuning voor dynamisch gedeeld geheugen aangeboden in de vorm van functies. Het zou voor de gebruiker intu¨ıtiever zijn indien een type was aangemaakt dat indexering door middel van de [] operator ondersteunt. Het probleem bij het definiëren van een nieuw type is echter dat, aangezien Julia dynamisch getypeerd is, dit leidt tot allerhande extra overhead, zoals encapsulatie van variabelen. Dit correct ondersteunen was te complex en tijdrovend om in het kader van deze masterproef te doen. Daarom is er voor een eenvoudiger alternatief met functies gekozen, dat nog steeds dezelfde functionaliteit aanbiedt. Deze functies passen, net zoals bij de intrinsics, rechtstreeks inlinen van LLVM instructies toe voor het definiëren en indexeren van gedeeld geheugen. Voor de implementatie van cuSharedMem() voor het definiëren van gedeeld geheugen schoot de implementatie van de inlining functionaliteit echter tekort. Instructies kunnen

28


immers enkel in functiebereik toegevoegd worden. Definitie van gedeeld geheugen vereist het toevoegen van een globale variabele in de LLVM module, wat vanuit functiebereik niet mogelijk is. Bijgevolg moest de inlining functionaliteit uitgebreid worden met een variant die rechtstreeks code toevoegt in de LLVM module, buiten functiebereik. De cuSharedMem() functie zal dus een globale variabele in de LLVM module toevoegen en vervolgens in de functie de nodige LLVM instructies inlinen om een referentie naar die globale variabele te verkrijgen. Om de verkregen referentie naar gedeeld geheugen vervolgens in een Julia variabele op te slaan, moeten ook nog extra conversie instructies gegenereerd worden. De addrspace annotaties die moeten toegevoegd worden zijn in LLVM een onderdeel van het datatype van een variabele. Zo is bijvoorbeeld het datatype van een gedeelde integer -rij int addrspace(3). Aangezien Julia geen rekening houdt met de verschillende adresruimtes op een GPU en enkel de generieke LLVM adresruimte ondersteunt, is het niet mogelijk om een datatype met addrspace annotatie in een Julia variabele op te slaan. De NVTPX backend voorziet echter intrinsics om geheugenreferenties te converteren van GPU adresruimtes naar de generieke adresruimte en terug. De functie cuSharedMem() converteert dus intern de referentie naar het aangemaakte gedeelde geheugen naar een generieke geheugenwijzer. De resulterende referentie kan dan wel in een Julia variabele opgeslagen worden. Deze variabele wordt dan aan de getCuSharedMem() en setCuSharedMem() functies doorgegeven, die intern dan de generieke geheugenwijzer omzetten naar een geheugenwijzer in de adresruimte voor gedeeld geheugen. Aangezien de implementatie van de getCuSharedMem() en setCuSharedMem() functies gebruikt maakt van het inlinen van LLVM instructies, wordt de toegang tot gedeeld geheugen niet door Julia codegeneratie verwerkt. Er kan dus geen gebruik gemaakt worden van de functionaliteit die de indexering startend vanaf 1 omzet in een indexering startend vanaf 0, zoals besproken in Paragraaf 3.3.4. Bijgevolg wordt van de opgegeven index van deze functies 1 afgetrokken, alvorens de waarde in de instructies te gebruiken. Om de allocatie van geheugen op de GPU en de transfer van data tussen CPU en GPU geheugen efficiënt te kunnen laten verlopen, is het nodig dat de gebruiker kan aanduiden welke kernel argumenten als input en/of output gebruikt worden. Hiervoor werden volgende wrappers voorzien, die gebruikt worden zoals beschreven in Paragraaf 3.1: CuIn: duidt aan dat een argument enkel inputwaarden voor de kernel bevat. Het is dus nodig om deze waarden naar het GPU geheugen te transporteren. Aangezien de programmeur dus geen nuttige resultaatwaarden in deze variabele verwacht, hoeft ze niet terug getransporteerd te worden. CuOut: duidt aan dat een argument de resultaatwaarden van de kernel zal bevatten.


29

Aangezien deze variabele dus initieel geen nuttige data bevat en door de kernel ingevuld zal worden, is het niet nodig data naar de GPU te sturen. Enkel de juiste hoeveelheid geheugen alloceren is voldoende. De resultaatwaarden moeten na afloop van de kernel natuurlijk wel terug getransporteerd worden. CuInOut: duidt aan dat een argument waarden bevat die nodig zijn voor de berekeningen in de kernel, maar dat deze waarden door de kernel gewijzigd zullen worden en nuttige eindresultaten zullen bevatten. Bijgevolg moet de data zowel naar de GPU als terug getransporteerd worden. Dit is nuttig bij bijvoorbeeld inplace berekeningen.

De functionaliteit implementeren die een indexering buiten de grenzen van een rij onderschept, vereist complexe wijzigingen. De gebruiker zal immers noch in de kernelargumenten, noch in de code van de kernel hiermee rekening houden. De kernel beschikt niet over de nodige informatie om de grenzen te kunnen controleren. Deze informatie ter beschikking stellen vereist dat de functiedeclaratie en kernelcode aangepast wordt zodat de benodigde informatie als kernel argument doorgegeven wordt en op de juiste manier in de kernel code verwerkt wordt. Er werd dus besloten om voorlopig enkel een eenvoudig geval te controleren, namelijk of de index minstens 1 is.

3.3.6

Linken met libdevice

Aangezien het niet mogelijk is om dynamisch functies in te laden op de GPU, is het ook niet mogelijk om de libdevice bibliotheek dynamisch te linken met de geschreven GPU code. Bijgevolg wordt deze bibliotheek statisch gelinked, wat ertoe leidt dat de volledige bibliotheek toegevoegd wordt aan de LLVM module die de GPU code bevat. Om te vermijden dat grote hoeveelheden code onnodig op de GPU worden ingeladen, moeten de functies die niet door de kernels gebruikt worden, na het linken terug uit de LLVM module gefilterd worden. De libdevice functies die door de kernels gebruikt worden moeten in de LLVM module gedeclareerd worden, net als andere functies die gebruikt worden. Voor het linken wordt een lijst opgesteld van alle aanwezige functies. Na het linken worden alle functies die niet in die lijst staan als intern gedeclareerd. Vervolgens wordt een GlobalOptimizerPass uitgevoerd die alle intern gedeclareerde functies verwijderd, met andere woorden, de overbodige libdevice functies. Om de prestaties te verhogen worden de gebruikte libdevice functies ge¨ınlined, teneinde extra functieoproepen te vermijden. Zoals vermeld in Paragraaf 2.1.3, beschikt een GPU ook over SFU’s voor het uitvoeren van een aantal complexe wiskundige functies, zoals sinus en cosinus. Gebruik makend van deze eenheden kunnen wiskundige functies sneller berekend worden, ten koste van de

30


accuraatheid. In grafische context is deze verhoogde prestatie belangrijker dan het verlies aan precisie, terwijl voor wetenschappelijke berekeningen en simulaties de accuraatheid primeert. Bijgevolg is gekozen om ondersteuning te bieden voor het gebruik van de libdevice bibliotheek.

3.3.7

Een in Julia geschreven kernel opstarten

Zoals te zien is in Paragraaf 3.1, is gekozen om het opstarten van een kernel via een macro te laten verlopen. Een macro in Julia kan gebruikt worden om code op arbitraire wijze te transformeren en expressies toe te voegen op de plaats in de syntaxis boom waar de macro voorkomt. Op deze manier kan de gebruiker voor het opstarten van een kernel wel de standaard syntaxis voor een functieoproep gebruiken. De @cuda macro vervangt deze dan door de correcte expressies voor het uitvoeren van volgende taken: Op basis van de wrappers, beschreven in Paragraaf 3.3.5, wordt voor de argumenten het nodige geheugen gealloceerd en de data getransporteerd. Aangezien de standaard syntaxis voor een functieoproep door de macro vervangen wordt, zal de Julia compiler deze kernel niet impliciet compileren naar LLVM IR. Bijgevolg moet deze compilatiestap expliciet door de macro getriggered worden. Vervolgens moet de LLVM IR omgezet worden naar PTX instructies die dan in de GPU geladen worden. Hoewel Julia de mogelijkheid biedt om van een functie de machinecode op te vragen, is dit onvoldoende: de GPU verwacht de PTX code van de volledige LLVM module die de kernels bevat, niet van één enkele kernel. Bijgevolg werd een functie toegevoegd aan Julia die de mogelijkheid biedt om in Julia code de PTX code van de voor de GPU bestemde LLVM module op te vragen. Vervolgens wordt aan de hand van de CUDA driver API deze PTX code in de GPU geladen en een functiereferentie naar de gewenste kernel opgevraagd. Deze referentie wordt gecached om het onnodig herhalen van deze compilatiestappen te vermijden. Gebruik makend van de CUDA driver API wordt de juist blok- en roosterconfiguratie ingesteld en wordt de kernel zelf opgestart. Nadat de kernel beëindigd is, moeten de resultaatwaarden terug naar het CPU geheugen getransporteerd worden. Opnieuw worden de kernelargumenten overlopen en wordt op basis van de wrapper het nodige datatransport uitgevoerd. Tenslotte wordt ook het geheugen dat door de macro op de GPU werd gealloceerd, vrijgegeven.

In de toekomst zouden deze taken ook door de Julia compiler uitgevoerd kunnen worden. De Julia compiler kan dan aangepast worden zodat deze taken uitgevoerd worden bij het


31

verwerken van een standaard functieoproep, wat betekent dat de gebruiker een kernel zou kunnen oproepen zoals een gewone functie. De programmeur moet bij deze oproep echter ook de gewenste blok- en roosterconfiguratie kunnen opgeven. Dit kan dan bijvoorbeeld aan de hand van de @cuda macro gebeuren. Een andere mogelijkheid is om de Julia syntaxis van een functieoproep uit te breiden met bijvoorbeeld de <<<...>>> syntaxis, zoals gedaan is voor CUDA C.

Hoofdstuk 4 Abstracties In dit hoofdstuk worden de abstracties voorgesteld die gebruikt zullen worden om de productiviteit van de programmeur te verhogen bij het gebruik van een GPU. Hiervoor worden eerst in Paragraaf 4.1 de verhouding tussen productiviteit en prestatie van de bestaande implementaties van de tracetransformatie geëvalueerd. Vervolgens worden de gekozen abstracties voorgesteld in Paragraaf 4.2 en wordt hun implementatie uitgelegd in Paragraaf 4.3. Tenslotte wordt in Paragraaf 4.4 verteld hoe deze abstracties concreet in de tracetransformatie werden toegepast.

4.1

Bestaande implementaties van de tracetransformatie

Algoritmes worden, afhankelijk van de context waarin ze ge¨ımplementeerd worden, door middel van verschillende talen of omgevingen ge¨ımplementeerd. Een typische taal voor veel wetenschappelijk onderzoek is MATLAB en zo is ook de initiële implementatie van de tracetransformatie in MATLAB gebeurd. Het aantal lijnen code nodig voor deze implementatie en de prestatie ervan is te zien in Tabel 4.1. Hoewel MATLAB na jaren ontwikkeling en optimalisatie zeer goede prestaties kan leveren, is dit voor veel implementaties niet het geval. Om de prestatie te verhogen kan, zoals bij veel hoogniveau programmeertalen, een tweestaps ontwikkelingsstrategie gehanteerd worden. Dit houdt in dat een algoritme in een hoogniveau taal geschreven wordt en de slecht presterende delen ervan opnieuw ge¨ımplementeerd worden in een laagniveau taal. In MATLAB is een tweestaps ontwikkelingsstrategie mogelijk door gebruik te maken van MEX-bestanden. Deze bevatten C, C++ of Fortran code en een toegangsfunctie mexFunction die door MATLAB opgeroepen wordt. Een voorbeeld voor het bereke32

4.2. KEUZE VAN ABSTRACTIES

33

nen van de som van twee vectoren aan de hand van een C-functie is te zien in Codevoorbeeld 4.1. In het geval van de tracetransformatie werd dit principe toegepast, wat inderdaad een prestatieverhoging met zich meebracht (zie Tabel 4.1). Het probleem hierbij is dat de programmeur de C taal moet kennen en rekening moet houden met alle technische details waarmee hij geconfronteerd wordt. Daarnaast wordt de implementatie van het algoritme verspreid over twee talen, wat de leesbaarheid zeker niet ten goede komt. Tenslotte brengt deze tweestaps ontwikkelingsstrategie vaak overhead met zich mee, zoals te zien is in Codevoorbeeld 4.1. Om te onderzoeken welke mogelijkheden er zijn om zowel productiviteit als prestatie te bieden werd de tracetransformatie opnieuw ge¨ımplementeerd in Julia. Dit omdat Julia, in tegenstelling tot MATLAB, opensource is en dus toegang biedt tot de interne werking van de taal en het compilatieproces. Zo kan de taal eenvoudiger geanalyseerd, geoptimaliseerd en uitgebreid worden. Om de prestatie van de resulterende code te evalueren werd ook een versie van de tracetransformatie ge¨ımplementeerd die als basisprestatie geldt. Hiervoor werd gebruik gemaakt van C++ voor de volledige implementatie, zonder tweestaps ontwikkelingsstrategie. Daarenboven is C++ ook een zeer populaire taal waar veel “accelleratoren” en uitbreidingen voor worden aangeboden. Een voorbeeld daarvan is CUDA, waarmee NVIDIA GPU’s geprogrammeerd kunnen worden. Aangezien bepaalde stappen van de tracetransformatie, zoals het uitvoeren van een functionaal op alle kolommen, zich zeer goed lenen tot parallelle uitvoering, levert het gebruik van een GPU een extra prestatiewinst op, wat te zien is in Tabel 4.1. Bij al deze implementaties is duidelijk te zien dat steeds een afweging gemaakt wordt tussen productiviteit en prestaties. MATLAB en Julia bieden, als hoogniveau talen, vooral een hoge productiviteit en boeten hiervoor in zekere mate in aan prestaties. Betere prestaties behalen kan door het toepassen van een tweestaps ontwikkelingsstrategie, zoals bij de MATLAB implementatie gebeurde, of het overschakelen naar een laagniveau taal, zoals C++. Maar dit resulteert dan weer in een negatieve invloed op de productiviteit. De beste prestaties worden behaald bij het gebruik van een GPU, maar de implementatiecomplexiteit is hier dan ook het hoogst.

4.2

Keuze van abstracties

De bedoeling van deze masterproef is te bepalen in welke mate productiviteit en prestaties verenigd kunnen worden, specifiek toegepast op het gebruik van GPU’s. Zoals besproken in Paragraaf 2.2 kan de productiviteit van de programmeur verhoogd worden

34

HOOFDSTUK 4. ABSTRACTIES

#include "mex.h" // nodig voor interactie met MATLAB void vadd(const double *a, const double *b, double *c, int len ,→ ) { int i; for(i=0;i
4.2. KEUZE VAN ABSTRACTIES Omgeving MATLAB MATLAB & C Julia C++ C++ & CUDA

35

Uitvoeringstijd T1 functionaal (s) Implementatiekost (SLOCs) 0.547 355 0.344 844 0.667 712 0.302 1489 0.020 2616

Tabel 4.1: Prestatie en productiviteit van de verschillende implementaties van de tracetransformatie

door hem abstracties aan te bieden. Deze stellen de programmeur in staat zich bezig te houden met de logica in een algoritme, zonder geconfronteerd te worden met technische implementatiedetails. Er zijn natuurlijk zeer veel verschillende abstracties uit allerhande contexten die gebruikt kunnen worden om de productiviteit van de programmeur te verhogen. In de context van een GPU is parallellisme fundamenteel, dus werd de aandacht gericht op abstracties waarvan de implementatie zich leent tot parallelle uitvoer. Rekening houdend met de case van de tracetransformatie is gekozen voor een aantal bekende abstracties, geselecteerd uit [24]: map, reduce en scan.

4.2.1

Map

Een map voert een bepaalde elementaire functie uit op elk element uit een collectie, zoals te zien is in Figuur 4.1. Belangrijk hierbij is dat de uitvoering van de elementaire functie op een element geen bijwerkingen heeft op de uitvoering van deze functie op de andere elementen. Deze onafhankelijkheid leidt tot het “embarrassingly parallel” karakter van deze abstractie: de verwerking van de verschillende elementen kan zonder veel moeite volledig parallel uitgevoerd worden. Er dient op geen enkel moment tussen de draden gesynchroniseerd te worden, tenzij op het einde van de functie. In de context van de tracetransformatie komt een map bijvoorbeeld overeen met het uitvoeren van een functionaal op de kolommen van een geroteerde figuur. De kolommen zijn dan de elementen van de abstractie en de functionaal is de elementaire functie. De toepassing van een functionaal op een kolom gebeurt onafhankelijk van de andere kolommen.

36


Figuur 4.1: Grafische weergave van een parallelle map

4.2.2

Reduce

Het parallel uitvoeren van de berekeningen van een functionaal introduceert afhankelijkheden. Zo berekent bijvoorbeeld de Radon transformatie de som van alle elementen in een kolom. Om dit parallel uit te kunnen voeren kan dus geen map gebruikt worden. Een passende abstractie voor deze transformatie is de reduce. Hierbij wordt een combinatiefunctie gebruikt om alle elementen uit een collectie paarsgewijs te combineren tot één resultaatwaarde. De volgorde waarin de elementen gecombineerd worden ligt niet vast. Hetzelfde resultaat uitkomen voor elke mogelijke volgorde vereist bepaalde algebra¨ısche eigenschappen van de combinatiefunctie: indien de combinatiefunctie associatief is, kunnen de verschillende paarsgewijze combinaties van elementen o.a. als volgt herordend worden:

s = (((((((a0 ⊗ a1 ) ⊗ a2 ) ⊗ a3 ) ⊗ a4 ) ⊗ a5 ) ⊗ a6 ) ⊗ a7 ) s = (((a0 ⊗ a1 ) ⊗ (a2 ⊗ a3 )) ⊗ ((a4 ⊗ a5 ) ⊗ (a6 ⊗ a7 ))) Zoals te zien in Figuur 4.2, ontstaat een boom van bewerkingen waarbij de operaties op hetzelfde niveau onafhankelijk van elkaar en dus in parallel uitgevoerd kunnen worden. Daarentegen is elk niveau afhankelijk van de resultaten uit het vorige niveau. De niveaus moeten dus serieel uitgevoerd worden, wat synchronisatie tussen de draden vereist.

4.2. KEUZE VAN ABSTRACTIES

37

Figuur 4.2: Grafische weergave van een parallelle reduce

4.2.3

Scan

Net als bij de Radon functionaal, ontstaan afhankelijkheden tussen de berekeningen van de T1 functionaal wanneer deze parallel uitgevoerd wordt. In tegenstelling tot de Radon transformatie echter, zijn bepaalde tussenresultaten van deze berekeningen onderdeel van de eindoplossing. De herordening van bewerking die gedaan werd in Paragraaf 4.2.2, is hier dan niet toepasbaar. Bijgevolg kon dus geen reduce gebruikt worden en is een derde abstractie ge¨ımplementeerd: de scan. Net als bij de reduce wordt een combinatiefunctie uitgevoerd op alle elementen uit een collectie, maar nu met behoud van bepaalde tussenuitkomsten: voor elk element uit de input collectie wordt er een resultaat berekend dat de combinatie is van alle elementen tot en met dat element. Ook deze abstractie kan door middel van een boomstructuur geparallelliseerd worden, maar aangezien de bepaalde tussenuitkomsten onderdeel zijn van het eindresultaat, zal de structuur complexer zijn dan in het geval van een reduce. Dit is te zien in Figuur 4.3. Opnieuw zijn de bewerkingen die in een niveau uitgevoerd worden onafhankelijk van elkaar, terwijl een niveau wel de resultaten van het vorige niveau nodig heeft.

38


Figuur 4.3: Grafische weergave van een parallelle scan

4.3. IMPLEMENTATIE VAN DE ABSTRACTIES

4.3

39

Implementatie van de abstracties

Bibliotheek in Julia De voornoemde abstracties zouden als een bibliotheek in CUDA ge¨ımplementeerd kunnen worden en via CUDA.jl dan in Julia ter beschikking gesteld worden. Een andere mogelijkheid is de abstracties in Julia te implementeren, waarbij de infrastructuur beschreven in Hoofdstuk 3, de abstracties dan omzet naar machinecode voor de GPU. Het belangrijkste voordeel van de tweede methode is dat er globale analyses en optimalisaties uitgevoerd kunnen worden die in het geval van een statisch gecompileerde CUDA bibliotheek niet mogelijk zijn. De Julia compiler past namelijk JIT compilatie toe en zal dus een functie pas compileren wanneer deze voor het eerst opgeroepen wordt. Informatie die enkel tijdens de uitvoering ter beschikking is, kan gebruikt worden om de gegenereerde code te optimaliseren. Tijdens compilatie van een CUDA bibliotheek is dergelijk informatie niet aanwezig. Een tweede voordeel dat de implementatie van de abstracties in Julia biedt is dat de conceptuele werking van de abstracties losgekoppeld wordt van het CUDA raamwerk. Bijgevolg is het mogelijk ook ondersteuning in te bouwen voor andere infrastructuren, in de context van GPU’s bijvoorbeeld voor OpenCL, zonder de implementatie van de abstracties te hoeven wijzigen. Tenslotte zorgt het hoogniveau karakter van Julia ervoor dat het schrijven van code productiever verloopt dan in geval van CUDA, wat te zien is in Tabel 4.1. Dit geldt dus ook bij het implementeren van deze en toekomstige abstracties.

Concrete implementatie Voor de implementatie van de map kan sterk gesteund worden op het SIMT principe dat een GPU aanbiedt: meerdere draden voeren dezelfde instructies uit op verschillende data elementen. Dit komt overeen met hoe een map conceptueel werkt. Bijgevolg wordt per inputelement een draad voorzien. Aangezien de draden onafhankelijk van elkaar hun element verwerken, kunnen de draden verdeeld worden over meerdere blokken, teneinde de hardware resources op de GPU optimaal te benutten. De draaden blokconfiguratie wordt dan ook aangepast op basis van de grootte van de inputcollectie. De scan werd ge¨ımplementeerd volgens het principe dat te zien is in Figuur 4.3. Omwille van de data-afhankelijkheden die hierin aanwezig zijn, worden de draden die de scan uitreken best in hetzelfde blok geplaatst. Op die manier kan namelijk gedeeld geheugen gebruikt worden om tussenuitkomsten op te slaan, wat sneller is dat het globaal geheugen. Deze implementatie is niet de meest optimale implementatie. In [25] wordt een

40


efficiënter alternatief voorgesteld, dat echter complexer is om te implementeren. Initiële tests hiervan leverden slechtere prestatie op dan het suboptimale algoritme en de implementatie verder optimaliseren is een tijdrovende bezigheid. Aangezien het vinden van de beste implementatie voor de scan niet het doel van deze masterproef is, is gekozen om gebruik te maken van het eenvoudigere algoritme dat, hoewel niet optimaal, toch een significante prestatiewinst oplevert. Voor de implementatie van de reduce wordt hetzelfde algoritme als voor de scan gebruikt, maar wordt enkel de totaalwaarde teruggegeven, niet de tussenuitkomsten. Deze implementatie zal, omwille van de complexere boomstructuur, meer bewerkingen uitvoeren dan strikt noodzakelijk is voor de reduce. Dit is te zien wanneer Figuur 4.2 en Figuur 4.3 vergeleken worden. In de praktijk is deze meerkost echter beperkt, aangezien het aantal niveaus in de boomstructuur er niet door toeneemt. Dit is eveneens te zien in Figuur 4.2 en Figuur 4.3. Ook voor de reduce zijn er beter presterende, maar complexere alternatieven [26]. De keuze om deze niet te gebruiken is gebaseerd op dezelfde overwegingen als bij de implementatie van de scan. Rekening houdend met het feit dat deze abstracties op een GPU worden uitgevoerd, bevatten voornoemde implementatie van de scan en reduce een beperking die optreedt wanneer deze abstracties toegepast worden op veel verschillende collecties. Dit is bijvoorbeeld het geval bij de tracetransformatie, waar een scan of reduce toegepast wordt op elke kolom van een figuur. De huidige implementatie leidt ertoe dat voor elke kolom een kernel opgestart moet worden, wat een significante kost met zich meebrengt. Daarnaast beperkt de uitvoering van de reduce en scan zich, zoals hiervoor vermeld, tot slechts één blok, terwijl een GPU beschikt over hardware die om meerdere blokken tegelijk uit te voeren. Bijgevolg is gekozen om de implementatie van de reduce en scan als volgt te wijzigen. In plaats van slechts één collectie van elementen als argument op te geven, wordt een verzameling van collecties opgeven als een collectie waarvan de elementen subcollecties zijn. Daarnaast worden de reduce en de scan als elementaire functie opgenomen in een map. De map zal dan zijn elementaire functie op de elementen van de collectie uitvoeren, in dit geval dus de reduce of scan op de verschillende subcollecties. Op deze manier wordt de granulariteit van de oproepen en dus ook het aantal oproepen verlaagd. Voor de tracetransformatie wil dit zeggen dat de programmeur niet langer voor elke kolom een reduce of scan moet oproepen. Voor één figuur wordt nu één gewijzigde reduce of scan opgeroepen die als collectie de volledige figuur meekrijgt, waarbij de verschillende kolommen als subcollecties beschouwd worden. Dit principe wordt, voor het geval van een scan, schematisch weergegeven in Figuur 4.4.

4.3. IMPLEMENTATIE VAN DE ABSTRACTIES

Figuur 4.4: Grafische weergave van de gewijzigde scan

41

42

4.4


Gebruik van de abstracties in de tracetransformatie

De abstracties uit Paragraaf 4.2 kunnen als volgt gebruikt worden om de Radon, T1 en T2 functionalen uit de tracetransformatie te implementeren. De Radon functionaal berekent de som van de elementen in een kolom, en dit voor elke kolom in een figuur: # rijen in kolom

radon(kolom) =

X

kolom[i]

(4.1)

i=1

Deze berekening komt exact overeen met het uitvoeren van een reduce op een kolom met de optelling als combinatiefunctie. Aangezien deze reductie op elke kolom onafhankelijk dient berekend te worden, werd de aangepaste reduce toegepast. Bij de T1 functionaal wordt eveneens een sommatie uitgevoerd op elke kolom van de figuur, maar de berekening is complexer dan bij de Radon functionaal: # rijen in kolom

X

T 1(kolom) =

i=m

kolom[i] ∗ (i − m)

(4.2)

Hierbij is m de index van de gewogen mediaan van de waarden in de kolom. Deze mediaanberekening werd aan de hand van twee abstracties ge¨ımplementeerd. Eerst wordt een scan toegepast om de totaalwaarde en tussenuitkomsten van een kolom te kennen. Vervolgens wordt de index van de mediaan bepaald als i zodanig dat: gescande kolom[i − 1] <

totaal totaal ∧ gescande kolom[i] ≥ 2 2

Aangezien dit onafhankelijk bewerkingen zijn, zou een map gebruikt kunnen worden. Dit is volgens de zuivere definitie van een map echter niet correct omdat meer dan 1 inputelement nodig is voor het uitvoeren de elementaire functie. Hiervoor kan in principe een stencil gebruikt worden, maar aangezien dit slechts een veralgemening van de map is met meerdere inputelementen, werd dit dan ook als een map ge¨ımplementeerd. De berekening van Vergelijking (4.2) gebeurt eveneens in twee stappen. Eerst moeten alle elementen vermenigvuldigd worden met (i − m), waarvoor een map kan gebruikt worden. Het resultaat daarvan wordt gebruikt als input voor de reduce met de optelling als combinatiefunctie. Aangezien deze echter pas mag beginnen rekenen vanaf index m, zal voorgaande map de elementen met index i < m vermenigvuldigen met 0. De T2 functionaal is zeer gelijkaardig aan de T1 functionaal: # rijen in kolom

T 2(kolom) =

X i=m

kolom[i] ∗ (i − m)2

(4.3)

4.4. GEBRUIK VAN DE ABSTRACTIES IN DE TRACETRANSFORMATIE

43

Hierbij is m opnieuw de index van de gewogen mediaan van de waarden in de kolom. De berekening van de gewogen mediaan is volledig dezelfde als bij de T1 functionaal en gebruikt dus dezelfde abstracties. De berekening van Vergelijking (4.3) gebruikt eveneens dezelfde abstracties als de de berekening van Vergelijking (4.2). Alleen zal de map alle elementen nu vermenigvuldigen met (i − m)2 in plaats van (i − m).

Hoofdstuk 5 Resultaten en analyse Gebruik makend van de in Hoofdstuk 3 ontwikkelde infrastructuur en de abstracties uit Hoofdstuk 4 kan nu geëvalueerd worden in welke mate de prestatie van de tracetransformatie in Julia verbeterd kan worden door gebruik van de GPU, zonder in te boeten op vlak van productiviteit. Deze evaluatie verloopt in meerdere stappen waarbij het vertrekpunt de beste prestatie biedt, maar de laagste productiviteit. Elke volgende stap probeert de productiviteit te verhogen met een zo klein mogelijke impact op de prestaties. In Paragraaf 5.2 wordt voor elke stap de invloed op de prestaties geëvalueerd. Daarna wordt in Paragraaf 5.3 gekeken wat de eigenlijke productiviteitswinst van elke stap is. Tenslotte wordt in Paragraaf 5.4 een aantal interessante mogelijkheden voor verder onderzoek vermeld.

5.1

Informatie over de opstelling

De experimenten gebeurden op basis van de tracetransformatie, meer specifiek de Radon, T1 en T2 functionaal. De Radon functionaal kan, zoals te zien in Paragraaf 4.4, triviaal berekend worden aan de hand van een reduce en wordt dan ook niet verder behandeld in deze tekst. De T2 functionaal is, zoals opnieuw te zien in Paragraaf 4.4, wat betreft het gebruik van de gedefinieerde abstracties volledige gelijk aan de T1 en is dan ook in deze masterproef niet verder uitgewerkt omdat de T2 functionaal ten opzicht van de T1 functionaal geen meerwaarde biedt. Bijgevolg hebben alle resultaten in dit hoofdstuk enkel betrekking op de T1 functionaal. Tenzij waar anders vermeld, is de figuur waarop de tracetransformatie wordt toegepast afkomstig uit [27]. Het is een realistische figuur in zwart-wit van 150 × 150 pixels, getrokken door een bewakingscamera in een tunnel [28]. Deze figuur is te zien in Figuur 5.1. De delen van de tracetransformatie die op de GPU uitgevoerd worden, zijn de rotatie van 44

5.2. PRESTATIEANALYSE

45

Figuur 5.1: Figuur waarop de tracetransformatie werd toegepast

de figuur en de berekening van de T1 functionaal. Voor alle meetwaarden werd de tracetransformatie steeds 30 keer uitgevoerd. Elke iteratie van de tracetransformatie voert de T1 functionaal uit op elke rotatie van de figuur gaande van 1° tot 359° in stappen van 1°. Hoewel de infrastructuur uit Hoofdstuk 3 automatisch beheer van het GPU geheugen aanbiedt, gebeurt dit momenteel op na¨ıeve beslissingen. Omdat dit een grote invloed heeft op de totale uitvoeringstijd van de tracetransformatie, werd voor die metingen het geheugen manueel beheerd. Meer informatie hierover is te vinden in Paragraaf 5.4. De hardware waarop de experimenten werden uitgevoerd omvat een Intel Core i7-3770K (Quad core @ 3.5GHz) CPU en een NVIDIA GeForce GTX Titan GPU, het besturingssysteem dat daarop draait is Debian unstable (“jessie”) op basis van de Linux kernel versie 3.13-1-amd64 x86 64. Aangezien Julia nog actief in ontwikkeling is, werd voor het bouwen van de infrastructuur en het uitvoeren van de metingen gebruik gemaakt van de op dat moment meest recente versie, namelijk 0.3.0-prerelease+3033 (2014-05-17 15:19 UTC). De gebruikte versie van CUDA is versie 5.5 en de gebruikte versie van LLVM is versie 3.3.

5.2

Prestatieanalyse

Figuur 5.2 geeft een algemeen beeld van de prestatie van de verschillende implementatiestappen. In Tabel 5.1 worden de exacte meetwaarden vermeld, terwijl de gemiddelde meetwaarden en relatieve versnelling in Tabel 5.2 worden weergegeven. De beste prestatie werd behaald met de implementatie van de tracetransformatie in C++ & CUDA en deze wordt dus als referentiepunt gebruikt. De andere implementaties passen verschillende methodes toe om in Julia gebruik te maken van de GPU. In wat volgt wordt de prestatie van deze methodes in detail geëvalueerd.

46

HOOFDSTUK 5. RESULTATEN EN ANALYSE

0.71 0.69 0.67

Uitvoeringstijd (s)

0.65 0.032 0.028 0.024 0.02

Julia

A CUD

A PTX CUD → & a i l a i Ju Jul

es

racti Abst

Figuur 5.2: Uitvoeringstijd van de tracetransformatie met T1 functionaal voor de verschillende implementaties

Omgeving Julia C++ & CUDA Julia & CUDA Julia → PTX Abstracties

Uitvoeringstijd in seconden Minimum Eerste kwartiel Mediaan Tweede kwartiel Maximum 0.6520 0.6592 0.6717 0.6752 0.7094 0.0190 0.0191 0.0191 0.0206 0.0228 0.0192 0.0193 0.0194 0.0208 0.0230 0.0240 0.0241 0.0242 0.0247 0.0288 0.0295 0.0297 0.0299 0.0301 0.0329

Tabel 5.1: Exacte meetresultaten van de uitvoeringstijd van de tracetransformatie met T1 functionaal

Omgeving Julia C++ & CUDA Julia & CUDA Julia → PTX Abstracties

Gemiddelde uitvoeringstijd (s) 0.6668 0.0201 0.0203 0.0248 0.0302

Gemiddelde relatieve versnelling 1.0 33.17 32.85 26.89 22.08

Tabel 5.2: Gemiddelde uitvoeringstijd van de tracetransformatie met T1 functionaal en gemiddelde relatieve versnelling t.o.v. de implementatie in Julia


5.2.1

47

Manueel gebruik van CUDA in Julia

Een eerste mogelijkheid om in Julia gebruik te maken van de GPU, is door toepassing van een tweestaps ontwikkelingsstrategie. Hierbij wordt, zoals vermeld in Hoofdstuk 3, gebruik gemaakt van het pakket CUDA.jl. Om zo goed mogelijk de meerkost van deze tweestaps ontwikkelingsstrategie te evalueren, wordt voor de verschillende kernels dezelfde code gebruikt als voor de referentie implementatie in C++ & CUDA. Zoals te zien is in Figuur 5.2, lijkt er op het eerste zicht geen verlies van prestatie op te treden. Na de uitvoering van de code grondiger onderzocht te hebben, kwam er wel een bron van vertraging aan het licht, namelijk bij het opstarten van een kernel. Zowel de CUDA runtime API, wat gebruikt wordt in CUDA C, als CUDA.jl bieden hiervoor semantisch gelijkaardige functies aan en maken hierin gebruik van de CUDA driver API. Het opstarten van kernels duurt bij CUDA C gemiddeld 3.46 × 10−6 seconden, terwijl dit bij CUDA.jl gemiddeld 5.51×10−6 seconden en dus bijna 60% trager is. De reden hiervoor is dat de argumenten voor de kernel naar de CUDA driver API doorgespeeld moeten worden als één enkel void** type. Intern zullen de opstartfuncties de argumenten dus samennemen tot een void**, wat in het geval van Julia encapsulatie van variabelen vereist en trager is. De reden dat dit uiteindelijk niet leidt tot een vertraging van de totale uitvoeringstijd van de tracetransformatie is de asynchrone uitvoer van kernels op de GPU. Nadat de kernel is opgestart, kan de CPU immers verder werken, terwijl de GPU de kernel uitvoert. De CPU kan dan nog steeds nieuwe kernels opstarten, welke dan door de GPU uitgevoerd worden van zodra de vorige afgewerkt is. Bij de tracetransformatie en ook andere algoritmes, is het echter zo dat de GPU zeer veel kernels uitvoert, terwijl de CPU niet veel meer doet dan deze kernels opstarten. De CPU slaagt er bij de C++ implementatie in om kernels naar de GPU te sturen voordat de vorige kernels berekend zijn en dus de GPU continu van werk te voorzien. De CPU zal zijn werk voltooien voordat de GPU klaar is en de totale uitvoeringstijd van het algoritme wordt dus bepaald door de uitvoeringstijd op de GPU. Bij de Julia implementatie, die gebruik maakt van CUDA.jl, is er, zoals hiervoor vermeld, meer vertraging tussen het opstarten van de verschillende kernels. Ondanks deze vertraging, slaagt de CPU er echter nog steeds in om de kernels snel genoeg op te starten en zijn taak sneller te voltooien dan de GPU. De CPU code in geval van de Julia implementatie is dus trager dan die van de C++ implementatie, maar dit verschil wordt gemaskeerd door de rekentijd op de GPU. Naarmate de te verwerken figuur kleiner is, heeft de GPU minder werk en zal het verschil in uitvoeringstijd op de CPU wel merkbaar worden, zoals te zien is in Figuur 5.3. De prestatie voor beide implementaties is gelijk bij een voldoende grote figuur, bij een resolutie van 64 × 64 is echter wel een verschil te zien.


Uitvoeringstijd (s)

48

10−1

Abstracties Julia → PTX Julia & CUDA C++ & CUDA

10−2 64x64 128x128 256x256 512x512 Resolutie van de figuren (pixels) Figuur 5.3: Gemiddelde uitvoeringstijd van de T1 functionaal voor verschillende resoluties van een realistische afbeelding voor gezichtsherkenning afkomstig van een bewakingscamera

5.2.2

Compilatie van Julia naar PTX

Een alternatieve methode die de tweestaps ontwikkelingsstrategie vermijdt, is het gebruiken van de in dit onderzoek ontwikkelde infrastructuur die het mogelijk maakt om Julia code automatisch te compileren naar PTX instructies, zoals beschreven in Hoofdstuk 3. De resulterende uitvoeringstijd is opnieuw te zien in Figuur 5.2. Om te bepalen wat de vertraging van de code veroorzaakt, dienen twee aspecten onderzocht te worden. Ten eerste moet bepaald worden hoe goed de Julia code naar PTX instructies gecompileerd wordt. Daarnaast dient, net als in Paragraaf 5.2.1, bekeken te worden in welke mate de CPU erin slaagt om tijdig nieuwe taken naar de GPU te sturen. Om de kwaliteit van de compilatie van Julia code naar PTX instructies te onderzoeken, werd gekeken naar de uitvoeringstijd van de resulterende kernels op de GPU. De relatieve uitvoeringstijden ten opzichte van de CUDA C implementatie is te zien in Tabel 5.3. Dit toont aan dat de compilatie redelijk goed verloopt, maar dat er toch nog ruimte voor verbetering is. Een duidelijk prestatieverschil is te zien bij de rotate kernel en de find weighted median kernel. De belangrijkste reden hiervoor is het prestatieverschil tussen 64bit en 32bit vlottende komma berekeningen dat inherent is aan de GPU. In de rotate en find weighted median kernels wordt immers een deling door een constante uitgevoerd. Julia bepaalt het datatype van deze constante op basis van de architectuur van de CPU, in dit geval 64bit. Dit resulteert in een 64bit deling die dus op de GPU slechter presteert dan een 32bit deling. Wanneer expliciet wordt opgegeven dat het bij deze constante gaat om een 32bit vlottende komma getal, wordt de 64bit deling vermeden en is de relatieve


49

Kernel Relatieve uitvoeringstijd rotate 1.52 prescan 0.98 find weighted median 1.74 1.08 t1 Tabel 5.3: Relatieve uitvoeringstijd van de Julia kernels ten opzichte van de CUDA C kernels

uitvoeringstijd voor de rotate functie nog slechts 1.17 en voor de find weighted median functie nog slechts 1.13. De programmeur dit expliciet laten aanduiden is echter niet gewenst, aangezien hij dan geconfronteerd wordt met hardware eigenschappen. Dit is juist wat men in een productieve hoogniveau taal zoveel mogelijk probeert te vermijden. De keuze tussen 32bit en 64bit moet dus door de Julia compiler gemaakt worden. Om conform te zijn met de huidige implementatie en documentatie van Julia, zou dit moeten gebeuren op basis de architectuur van de GPU. Voornoemd probleem is daarmee echter niet opgelost, aangezien bijna alle moderne GPU’s een 64bit architectuur hebben. Een andere oplossing is om de Julia compiler zoveel mogelijk voor 32bit bewerkingen te laten kiezen. Belangrijk is dan dat de gewenste precisie van de berekeningen niet in het gedrang komt en 64bit getallen ten onrechte omgezet worden naar 32bit. Om Julia code, die bedoeld is voor de GPU, ook uit te voeren op de GPU, dient gebruikt gemaakt te worden van de @cuda macro zoals beschreven in Hoofdstuk 3. Deze macro voert extra administratieve taken uit, die een vertraging met zich meebrengen. Zo moet bepaald worden voor welke functieargumenten geheugen op de GPU dient gealloceerd te worden en welke datatransfers dienen uitgevoerd te worden. Daarnaast worden de gecompileerde functies ook gecached en moet bij het oproepen die functie uit de cache opgehaald worden. Tenslotte start de macro de kernel op, net als in Paragraaf 5.2.1, wat hier ook dezelfde vertragingen introduceert. In totaal duurt het opstarten van de kernel aan de hand van de @cuda macro nu 1.58 × 10−5 seconden. De vertraging is nu, in tegenstelling tot bij Paragraaf 5.2.1, te groot zodat de GPU op sommige momenten geen werk heeft. Dit kan nagegaan worden met de NVIDIA profiler, die de uitvoering van kernels op een GPU analyseert. In de tijdlijn van de GPU, weergegeven in Figuur 5.4(a), is te zien dat de kernels elkaar onmiddellijk opvolgen; de korte tussenruimte is te wijten aan het configureren van de nieuwe kernel. In Figuur 5.4(b) daarentegen volgen de kernels elkaar niet snel genoeg op en wordt de GPU dus niet maximaal belast. Opnieuw is deze invloed te zien in Figuur 5.3. Wanneer de resolutie te klein is (64 × 64 of 128×128) kan de CPU de GPU niet bijhouden. De uitvoeringstijd wordt dan bepaald door

50


(a) Julia & CUDA - kernels volgen elkaar zo snel mogelijk op en de GPU wordt maximaal belast

(b) Julia → PTX - kernels volgen elkaar niet altijd snel genoeg op waardoor de GPU soms op werk moet wachten

Figuur 5.4: Tijdlijnen van de uitvoering van kernels op de GPU

de CPU en is onafhankelijk van de resolutie van de figuur. Voor te kleine afbeeldingen wordt dus een ondergrens voor de prestatie bereikt. Voor grotere figuren zal de GPU opnieuw meer werk hebben dan de CPU en wordt de vertraging ten opzichte van de vorige implementatie verklaard door het prestatieverschil van de kernels zelf (zie Tabel 5.3).

5.2.3

Gebruik van abstracties

Een laatste methode om de GPU in Julia te introduceren is door gebruik te maken van abstracties. De toepassing van de abstracties is gebeurd zoals beschreven in Paragraaf 4.4. Aangezien de rotatie functie ook een groot aandeel in het rekenwerk heeft, wordt deze eveneens op de GPU uitgevoerd. Door de complexe afhankelijkheden tussen input en output elementen is het echter niet mogelijk om deze te implementeren op basis van de abstracties uit Hoofdstuk 4. Dat niet alle berekeningen altijd in abstracties kunnen omgezet worden zal in de praktijk zeker voorkomen, dus is er gekozen om de rotatie functie ook hier manueel te implementeren, gebruik makend van dezelfde implementatie als in Paragraaf 5.2.2. Wat de productiviteit van de programmeur betreft is dit natuurlijk nadelig, maar dankzij de infrastructuur van Hoofdstuk 3 moet niet direct teruggevallen worden op CUDA C en wordt de tweestaps ontwikkelingsstrategie vermeden. Het gebruik van abstracties in de tracetransformatie leidt tot een vertraging in de uitvoeringstijd, zoals te zien is in Figuur 5.2. De implementatie van de abstracties maakt, net als in Paragraaf 5.2.2 gebruik van de @cuda macro, wat ertoe leidt dat taken opnieuw niet snel genoeg naar de GPU gestuurd kunnen worden. Ook hier is er dus voor kleinere figuren een ondergrens op de te behalen prestatie. De reden waarom het gebruik van abstracties echter nog trager is dan de implementatie van Paragraaf 5.2.2, is het feit dat de T1 functionaal bij gebruik van abstracties een extra

5.3. PRODUCTIVITEITSANALYSE Stap in berekening T1 Gewogen mediaan berekenen Sommatie

51

Julia functie gebruikte abstractie prescan scan find weighted median map map T1 reduce

Tabel 5.4: Functies en abstracties nodig voor het berekenen van de T1 functionaal bij manuele implementatie (Julia → PTX) en gebruik van abstracties

stap heeft ten opzichte van voornoemde implementatie. Dit is te zien in Tabel 5.4, waar aangegeven wordt welke functies of abstracties gebruikt worden voor welke berekeningen. Gebruik makende van abstracties wordt er dus een extra GPU kernel opgestart. Aangezien de prestatie beperkt wordt door de snelheid waarmee de kernels op de GPU worden opgestart, zorgt dit voor 25% extra vertraging. Een oplossing hiervoor is om de twee abstracties van de sommatie uit Tabel 5.4 tot combineren tot één enkele GPU kernel. Dit is zeker mogelijk, aangezien ze perfect serieel worden uitgevoerd: de output van de map is de input van de reduce en de combinatie van deze twee leidt dus tot de bekende mapreduce. Een nadeel van het voorzien van mapreduce is dat er op die manier zeer veel mogelijke combinaties van abstracties voorzien zouden moeten worden. Het zou beter zijn gebruik te maken van globale analyse om te detecteren dat deze abstracties gecombineerd kunnen worden. Dankzij de JIT compilatie in Julia kan de gegenereerde code dan hieraan aangepast worden.

5.3

Productiviteitsanalyse

In Paragraaf 5.2 werd gebruik gemaakt van de GPU om de prestatie van de Julia implementatie van de tracetransformatie te verhogen. Hierbij was het de bedoeling om de productiviteit van de programmeur zoveel mogelijk in stand te houden. Een mogelijke indicator hiervan is het aantal lijnen code dat nodig is om een algoritme te implementeren. Immers, wanneer dezelfde logica met minder lijnen ge¨ımplementeerd kan worden, zijn deze lijnen semantisch rijker en kan de programmeur op productievere wijze zijn algoritme omzetten in code. Het aantal lijnen code nodig voor de tracetransformatie en voor de Radon en T1 functionaal is te zien in Tabel 5.5. De eerste stap bij de integratie van de GPU in Julia, is het toepassen van een tweestaps ontwikkelingsstrategie, wat geëvalueerd werd in Paragraaf 5.2.1. Aangezien de logica van het algoritme nu in Julia geschreven wordt, gaat de productiviteit erop vooruit ten opzichte van de C++ implementatie. Deze logica vormt een groot deel van de totale implementatie, wat ertoe leidt dat er een grote reductie is in het totaal aantal lijnen code.

52

HOOFDSTUK 5. RESULTATEN EN ANALYSE Omgeving Julia C++ & CUDA Julia & CUDA Julia → PTX Abstracties

Volledige implementatie (SLOC) 604 1936 795 685 522

Radon (SLOC) 7 69 83 43 5

T1 (SLOC) 19 147 172 96 20

Tabel 5.5: Aantal lijnen code (SLOC) nodig voor het implementeren van de tracetransformatie, Radon functionaal en T1 functionaal, voor de verschillende omgevingen

Vergelijken we de productiviteit echter met de zuivere Julia implementatie, dan heeft de gehanteerde tweestaps ontwikkelingsstrategie een grote impact op de productiviteit. Dit resulteert in de grote toename in aantal lijnen code voor het implementeren van de Radon en T1 functionalen. Deze implementatie gebeurt nu immers in CUDA in plaats van in Julia. Ten opzicht van de implementatie in C++ & CUDA gaat de productiviteit erop vooruit, maar de tweestaps ontwikkelingsstrategie en het gebruik van CUDA zorgen voor een groot verschil in productiviteit ten opzichte van de implementatie in Julia.

Om de productiviteit te verhogen, werd in de volgende stap gebruik gemaakt van de infrastructuur van Hoofdstuk 3. Op deze manier kan de programmeur de kernels in Julia schrijven en wordt de tweestaps ontwikkelingsstrategie en de nood aan kennis van CUDA C vermeden. Dit heeft direct een grote invloed op de productiviteit van de programmeur. Wat het aantal lijnen code betreft is de implementatiekost voor de Radon en T1 functionaal gedaald ten opzichte van de de Julia & CUDA implementatie. Een vergelijkbare productiviteit als bij de implementatie in Julia zonder GPU is echter nog niet behaald.

In de laatste stap om de productiviteit te verhogen, wordt gebruik gemaakt van de abstracties, zoals beschreven in Paragraaf 4.4. De programmeur hoeft zich niet bezig te houden met de specifieke implementatie van de abstracties. Bepaalde technische aspecten van de GPU waar de programmeur in de vorige stap nog mee geconfronteerd werd, zoals het beheer van het gedeeld geheugen, blijven nu volledige verborgen voor de programmeur. De abstracties laten dus toe om op quasi transparante wijze gebruik te maken van een GPU. Zodoende kan de programmeur zich concentreren op de conceptuele werking van het algoritme en verloopt de implementatie zeer vlot. Dit is ook te zien in Tabel 5.5: gebruik makend van abstracties wordt dezelfde productiviteit als van de zuivere Julia implementatie behaald, maar dit met een significatie prestatiewinst.

5.4. TOEKOMSTIG WERK

5.4

53

Toekomstig werk

Hoewel de resultaten in Paragraaf 5.2 en Paragraaf 5.3 aantonen dat een combinatie van prestatie en productiviteit haalbaar is, zijn er op beide vlakken zeker nog mogelijkheden voor optimalisaties. Bij het gebruik van een GPU is er nog een factor die een significante invloed kan hebben op de prestatie, namelijk het transport van data tussen CPU en GPU geheugen. De @cuda macro die gebruikt wordt in Paragraaf 5.2.2 en Paragraaf 5.2.3 automatiseert het geheugenbeheer en datatransport is daar een onderdeel van. Bepalen wanneer data getransporteerd moet worden, gebeurt momenteel echter op na¨ıeve wijze. Inputdata voor een kernel zal steeds naar het GPU geheugen getransporteerd worden en outputdata ook steeds terug naar het CPU geheugen. Indien de inputdata nog niet aanwezig was op de GPU en de outputdata nodig is de CPU, zijn deze datatransfers nodig en wordt dus een correcte beslissing gemaakt. Het komt echter frequent voor dat meerdere kernels samenwerken voor een bepaalde rekentaak en dat de resultaten van de ene kernel enkel nodig zijn als input voor een andere kernel. Transfer van deze data naar het CPU geheugen is dan overbodig en zou vermeden moeten worden. De @cuda macro beschikt momenteel echter niet over de nodige globale informatie om dergelijke gevallen te herkennen. Alle tussenresultaten worden dus tussen CPU en GPU geheugen uitgewisseld. Deze situatie treedt ook op bij de tracetransformatie. Zo wordt bijvoorbeeld de rotatie van een figuur op de GPU uitgevoerd en is de geroteerde figuur enkel nodig voor de T functionalen die ook op de GPU worden uitgevoerd. Transfer van deze geroteerde figuur naar het CPU geheugen is dus overbodig, maar wordt door de huidige implementatie van de @cuda macro niet weggefilterd. Om te vermijden dat inefficiënt geheugenbeheer de prestaties van de infrastructuur uit Hoofdstuk 3 teveel be¨ınvloedt, werd in de @cuda macro nog de mogelijkheid voorzien om manueel het geheugen op de GPU te beheren. Bij het gebruik van abstracties moet de programmeur kunnen opgeven welke concrete elementaire of combinatiefunctie uitgevoerd moet worden. Momenteel wordt de code van de gewenste functies rechtstreeks in de abstractie geschreven en geeft de programmeur zijn keuze aan op basis van een extra argument van de abstractie. Hoewel dit een mogelijkheid is om verschillende functies aan de programmeur aan te bieden, beperkt het de flexibiliteit van de abstracties. Een beter alternatief dat de gebruiker volledige vrijheid geeft voor het beschrijven van de functie, is het gebruik van een expressie, bv ex = :(a+b). Deze expressie wordt als argument aan de abstractie doorgegeven en rechtstreeks in de abstractie ge¨ınlined. Hoewel dit weinig overhead geeft, worden de symbolen van de variabelen in de expressie rechtstreeks in de code toegevoegd. Het is dus belangrijk dat er voor die symbolen een conventie afgesproken wordt. Anderzijds wil dit ook zeggen dat er bij het gebruik van extra variabelen voor tussenbewerkingen naamconflicten kun-

54


nen optreden. Deze problemen kunnen eenvoudig vermeden worden door de bewerking aan de abstractie door te geven als een anonieme functie, bv (a,b)->a+b. Voor deze functie kan dan aparte LLVM IR gegenereerd worden, die daarna ge¨ınlined wordt door de compiler. Op die manier hoeven er geen naamconventies bepaald te worden en treden er ook geen naamconflicten op. Ook kunnen de experimenten die hier gebeurd zijn herhaald worden met andere toepassingen uit diverse wetenschapsdomeinen. Hiermee kan een sterkere bevestiging van de generaliteit van de gekozen abstracties bekomen worden. Uiteraard kunnen ook nieuwe abstracties onderzocht worden.

Hoofdstuk 6 Conclusie Het gebruik van een GPU kan de uitvoering van rekenalgoritmes in veel wetenschappelijk onderzoek sterk versnellen. Het gebruik van een GPU is echter complex en vereist kennis van technische en architecturale details die vaak ontbreekt. Is de kennis wel aanwezig, dan verloopt het integreren van een GPU in het algoritme nog vaak moeilijk omdat een tweestaps strategie gebruikt wordt: een deel van het algoritme wordt in een hoogniveau programmeertaal geschreven (b.v. Julia) terwijl de GPU-specifieke delen in een andere, laagniveau taal ontwikkeld worden (CUDA C). Dit alles bemoeilijkt vaak het gebruik van een GPU in wetenschappelijk onderzoek. In deze masterproef werd onderzocht hoe de programmeur op productieve wijze een GPU kan gebruiken om de prestatie van code te verhogen. Als oplossing hiervoor werden een aantal generieke abstracties gedefinieerd en ge¨ımplementeerd, en welke transparant op GPU uitgevoerd worden. De implementatie van de abstracties is gebeurd in Julia om de mogelijkheid te bieden de gegenereerde GPU code aan te passen op basis van informatie die enkel tijdens uitvoering ter beschikking is. Aangezien compilatie van Julia code naar GPU instructies niet ondersteund is in de huidige Julia compiler, werd deze uitgebreid met een infrastructuur die deze compilatie mogelijk maakt. Ook worden een aantal administratieve taken die gepaard gaan met het gebruik van een GPU door deze infrastructuur uitgevoerd, opnieuw om dit voor de programmeur te verbergen. Uit de resultaten van Hoofdstuk 5 is gebleken dat de prestatie van de tracetransformatie waarin gebruik gemaakt wordt van deze oplossing ongeveer 50% trager is dan het manueel programmeren en integreren van een GPU in het algoritme. Desalniettemin is de productiviteit van de programmeur bij gebruik van deze oplossing vele malen hoger dan bij de manuele integratie van een GPU. Bovendien is de prestatie nog steeds 22 keer hoger dan de gewone Julia implementatie waar geen gebruik gemaakt wordt van een GPU. Deze resultaten tonen aan dat de gekozen oplossing veelbelovend is. De bereikte implementatie ervan is echter zeker nog niet compleet, maar kan gezien worden als een 55

56

HOOFDSTUK 6. CONCLUSIE

prototype dat aantoont wat mogelijk is. Ook kan het gebruikt worden voor verder onderzoek naar andere abstracties, of naar manieren om de implementatie van de abstracties te verbeteren, b.v. door gebruikt te maken van globale optimalisatie.

Bibliografie [1] A. Remenyi, S. Szenasi, I. Bandi, Z. Vamossy, G. Valcz, P. Bogdanov, S. Sergyan, and M. Kozlovszky. Parallel biomedical image processing with gpgpus in cancer research. In 3rd IEEE International Symposium on Logistics and Industrial Informatics (LINDI), 2011, pages 245–248, Aug 2011. [2] Chia-Feng Juang, Teng-Chang Chen, and Wei-Yuan Cheng. Speedup of implementing fuzzy neural networks with high-dimensional inputs through parallel processing on graphic processing units. IEEE Transactions on Fuzzy Systems, 19(4):717–728, Aug 2011. [3] W. Vanderbauwhede and T. Takemi. An investigation into the feasibility and benefits of gpu/multicore acceleration of the weather research and forecasting model. In International Conference on High Performance Computing and Simulation (HPCS), 2013, pages 482–489, July 2013. [4] Jeff Bezanson, Stefan Karpinski, Viral B. Shah, and Alan Edelman. Julia: A fast dynamic language for technical computing. CoRR, abs/1209.5145, 2012. [5] E. Scott Larsen and David McAllister. Fast matrix multiplies using graphics hardware. In Proceedings of the 2001 ACM/IEEE Conference on Supercomputing, SC ’01, pages 55–55, New York, NY, USA, 2001. ACM. [6] Timothy J. Purcell, Ian Buck, William R. Mark, and Pat Hanrahan. Ray tracing on programmable graphics hardware. In Proceedings of the 29th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’02, pages 703–712, New York, NY, USA, 2002. ACM. [7] Chris J. Thompson, Sahngyun Hahn, and Mark Oskin. Using modern graphics architectures for general-purpose computing: A framework and analysis. In Proceedings of the 35th Annual ACM/IEEE International Symposium on Microarchitecture, MICRO 35, pages 306–317, Los Alamitos, CA, USA, 2002. IEEE Computer Society Press. 57

58

BIBLIOGRAFIE

[8] CUDA Nvidia. Programming guide, 2008. [9] Tianyi David Han and Tarek S Abdelrahman. hicuda: High-level gpgpu programming. IEEE Transactions on Parallel and Distributed Systems, 22(1):78–90, 2011. [10] Ross Ihaka and Robert Gentleman. R: A language for data analysis and graphics. Journal of Computational and Graphical Statistics, 5(3):299–314, 1996. [11] Eric Jones, Travis Oliphant, and Pearu Peterson. Scipy: Open source scientific tools for python. http://www. scipy. org/, 2001. [12] Stefan Van Der Walt, S Chris Colbert, and Gael Varoquaux. The numpy array: a structure for efficient numerical computation. Computing in Science & Engineering, 13(2):22–30, 2011. [13] Carl Friedrich Bolz, Antonio Cuni, Maciej Fijalkowski, and Armin Rigo. Tracing the meta-level: Pypy’s tracing jit compiler. In Proceedings of the 4th Workshop on the Implementation, Compilation, Optimization of Object-Oriented Languages and Programming Systems, ICOOOLPS ’09, pages 18–25, New York, NY, USA, 2009. ACM. [14] Floréal Morandat, Brandon Hill, Leo Osvald, and Jan Vitek. Evaluating the design of the r language. In James Noble, editor, ECOOP 2012 – Object-Oriented Programming, volume 7313 of Lecture Notes in Computer Science, pages 104–131. Springer Berlin Heidelberg, 2012. [15] Kurt Keutzer and Tim Mattson. A design pattern language for engineering (parallel) software. Intel Technology Journal, 13(4), 2009. [16] Krste Asanovic, Rastislav Bodik, James Demmel, Tony Keaveny, Kurt Keutzer, John Kubiatowicz, Nelson Morgan, David Patterson, Koushik Sen, John Wawrzynek, et al. A view of the parallel computing landscape. Communications of the ACM, 52(10):56–67, 2009. [17] Eric Holk, Milinda Pathirage, Arun Chauhan, Andrew Lumsdaine, and Nicholas D Matsakis. Gpu programming in rust: Implementing high-level abstractions in a systems-level language. In IEEE 27th International Parallel and Distributed Processing Symposium Workshops & PhD Forum (IPDPSW), 2013, pages 315–324. IEEE, 2013. [18] C. Lattner and V. Adve. Llvm: a compilation framework for lifelong program analysis transformation. In International Symposium on Code Generation and Optimization, 2004, pages 75–86, March 2004.

BIBLIOGRAFIE

59

[19] A. Kadyrov and M. Petrou. The trace transform and its applications. IEEE Transactions on Pattern Analysis and Machine Intelligence, 23(8):811–828, Aug 2001. [20] S. Srisuk, M. Petrou, W. Kurutach, and A. Kadyrov. A face authentication system using the trace transform. Pattern Analysis and Applications, 8(1-2):50–61, 2005. [21] A. Kadyrov and M. Petrou. Affine parameter estimation from the trace transform. IEEE Transactions on Pattern Analysis and Machine Intelligence, 28(10):1631–1645, Oct 2006. [22] SuhaibA. Fahmy, Christos-Savvas Bouganis, PeterY.K. Cheung, and Wayne Luk. Real-time hardware acceleration of the trace transform. Journal of Real-Time Image Processing, 2(4):235–248, 2007. [23] Lin Dahua. https://github.com/JuliaGPU/CUDA.jl, 2014. [24] Michael McCool, James Reinders, and Arch Robison. Structured parallel programming: patterns for efficient computation. Elsevier, 2012. [25] Hubert Nguyen. Gpu Gems 3. Addison-Wesley Professional, first edition, 2007. [26] Mark Harris et al. Optimizing parallel reduction in cuda. NVIDIA Developer Technology, 2:45, 2007. [27] Andres Frias-Velazquez, Carlos Ortiz, Aleksandra Pizurica, Wilfried Philips, and Gustavo Cerda. Object identification by using orthonormal circus functions from the trace transform. In 19th IEEE International Conference on Image Processing (ICIP), pages 2153–2156. IEEE, 2012. [28] Andrés Fr´ıas-Velázquez, Peter Van Hese, Aleksandra Piˇzurica, and Wilfried Philips. Vehicle classification for road tunnel surveillance. In Robert Paul Loce, Eli Saber, and Sreenath Rao Vantaram, editors, IS&T/SPIE Electronic Imaging, pages 86630M–86630M–6. International Society for Optics and Photonics, March 2013.

Lijst van figuren 2.1

Voorbeeld van een 3 × 4 blokconfiguratie en een 2 × 3 roosterconfiguratie

7

2.2

Organisatie van de CUDA runtime en de CUDA driver . . . . . . . . . .

9

2.3

Schematische weergave van de verschillende componenten op een GPU

.

10

2.4

Blokken van een kernel worden op dynamische wijze verdeeld over de aanwezig multiprocessors . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

11

2.5

Benchmark uitvoeringstijden relatief ten opzichte van C . . . . . . . . . .

14

2.6

Schematische weergave van de verschillende stappen in de tracetransformatie 17

4.1

Grafische weergave van een parallelle map . . . . . . . . . . . . . . . . .

36

4.2

Grafische weergave van een parallelle reduce . . . . . . . . . . . . . . . .

37

4.3

Grafische weergave van een parallelle scan . . . . . . . . . . . . . . . . .

38

4.4

Grafische weergave van de gewijzigde scan . . . . . . . . . . . . . . . . .

41

5.1

Figuur waarop de tracetransformatie werd toegepast . . . . . . . . . . . .

45

5.2

Uitvoeringstijd van de tracetransformatie met T1 functionaal voor de verschillende implementaties . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Gemiddelde uitvoeringstijd van de T1 functionaal voor verschillende resoluties van een realistische afbeelding voor gezichtsherkenning afkomstig van een bewakingscamera . . . . . . . . . . . . . . . . . . . . . . . . . .

48

Tijdlijnen van de uitvoering van kernels op de GPU . . . . . . . . . . . .

50

5.3

5.4

60

Lijst van tabellen 4.1

5.1 5.2 5.3 5.4 5.5

Prestatie en productiviteit van de verschillende implementaties van de tracetransformatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

Exacte meetresultaten van de uitvoeringstijd van de tracetransformatie met T1 functionaal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Gemiddelde uitvoeringstijd van de tracetransformatie met T1 functionaal en gemiddelde relatieve versnelling t.o.v. de implementatie in Julia . . .

46

Relatieve uitvoeringstijd van de Julia kernels ten opzichte van de CUDA C kernels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

Functies en abstracties nodig voor het berekenen van de T1 functionaal bij manuele implementatie (Julia → PTX) en gebruik van abstracties . . Aantal lijnen code (SLOC) nodig voor het implementeren van de tracetransformatie, Radon functionaal en T1 functionaal, voor de verschillende omgevingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

51

52

Parallelle abstracties voor het programmeren van GPU's in Julia

Recommend Documents