In modern systems CPUs as well as GPUs are equipped with multi-level memory architectures, where different levels of the hierarchy vary in latency and capacity. Therefore, various memory access models were studied. Such a model can be seen as an interface abstracting the user from the physical architecture details. In this paper we present a general and uniform GPU computation and memory access model based on bounded inhibitor Petri nets (PNs). Its effectiveness is demonstrated by comparing its throughputs to practical computational experiments performed with the usage of Nvidia GPU with CUDA architecture. Our PN model is consistent with the workflow of multithreaded GPU streaming multiprocessors. It models a selection and execution of instructions for each warp. The three types of instructions included in the model are: the arithmetic operation, the access to the shared memory and the access to the global memory. For a given algorithm the model allows to check how efficient the parallelization is, and whether a different organization of threads will improve performance. The accuracy of our model was tested with different kernels. As the preliminary experiments we used the matrix multiplication program and stability example created by Nvidia, and as the main experiment a binary version of the least significant digit radix sort algorithm. We created three implementations of the algorithm using CUDA architecture, differing in the usage of shared and global memory as well as organization of calculations. For each implementation the PN model was used and the results of experiments are presented in the work.
{anna.gogolinska,lukasz.mikulski,marcin.piatkowski}@mat.umk.pl
The inter-process communication over a common part of the memory shared by
processes is a usual performance bottleneck in multiprocessor environments. In
modern systems CPUs as well as GPUs are equipped with multi-level memory
architectures, where different levels of the hierarchy vary in latency and capacity.
By considering different local views of the processes on the common part of the
memory one can try to improve the processor utilization. Therefore, various
memory access models were studied, see for instance [
Another important issue is the task distribution between threads and CPU/GPU cores and the instruction scheduling, which can have a significant impact on the efficiency. Consider for instance running the three threads on a single processor depicted in Figure 1. Each of them performs a list of arithmetic operations interleaved by memory reads/writes consisting of a short preprocessing and then a longer period of waiting for the memory access (which is a usual situation in parallel computing). In the initial part of the computation each thread realizes its preprocessing for the memory access and then starts to wait for the access itself. This causes the processor idle period when all threads are waiting (marked as I). After that the arithmetic operations of all threads are executed simultaneously. They can be scheduled in such a way that one thread waits for the memory access while other threads perform their computations and there is no idle period (marked as II). Thanks to that the waiting period of a thread can be hidden behind the active computations of other threads.
t1: t2: t3: (I)
(II) – memory access – arithmetic operation
The main contribution of this paper is a general and uniform GPU
computation and memory access model based on bounded Petri nets [
The effectiveness of our model is demonstrated by comparing its
throughputs to practical computational experiments performed with the usage of Nvidia
GPU. We study the impact on the complexity of various parameters such as the
number of concurrent processes, and the level of memory used (shared memory
vs. global memory). The successful application of our model is discussed in
details, as a proof of concept, on a single example of digit radix sorting algorithm.
We do not want to discuss methods of parallelization, nor the most efficient way
of using different types of the memory. It is beyond the scope of this paper. Those
problems are very complex, and cannot be explained in such a short publication.
More about those topics can be found for instance in [
Our PN model is not the only one GPU efficiency model. There are quite
a few other GPU performance models, which can be divided into two groups
(according to [
(1) Calibrated performance models that make specific predictions and include
many lower levels of details. They usually contain many specialized parameters,
some of which may be difficult to obtain or calculate. The first example is a
model presented in [
(2) Asymptotic models for algorithm analysis at a high level of abstraction
that capture only the essential features of GPU architecture. One of the models
from this group is the model described in [
The paper is organized as follows. In the next section we describe Graphical Processing Units and CUDA Toolkit, focusing in particular on memory types (and organization). In Section 3 we recall some standard notions and notations related to Petri nets. Then we introduce our memory access and computation model followed by the description of radix sort algorithm. We also present the results of experiments conducted on the base of our implementations of the radix sort algorithm. We conclude the paper and give some directions of further research in the last section. 2
An intensive development of Graphical Processing Units (GPU in short) resulted
in construction of high performance computational devices, which besides the
graphical display management, allow also the execution of parallel general
purpose algorithms (not necessarily related to computer graphics). As opposed to
CPU consisting of a few cores optimized for sequential serial processing, GPU
contains thousands of smaller, more efficient cores designed for handling
multiple tasks simultaneously. The most popular ones are GPU’s produced by Nvidia
Corporation, supplied with Nvidia CUDA Toolkit (see [
A program running in heterogeneous environment equipped with GPU can be split into so-called host parts, which are executed by CPU, and so-called kernel parts, which are executed by GPU. The host part specify the kernel execution context and manages the data transfer between the host and the GPU memory. The kernel functions create a big number of threads allowing a highly parallel computation (where each thread runs the same kernel code). GPU threads, as opposed to CPU threads, are much lighter, hence their creation and controlling requires less CPU cycles.
All threads executed on GPU are organized into several equally-sized thread blocks, which in turn are organized into a grid.
A thread block is the set of concurrently executed threads. Such an execution and the thread cooperation can be coordinated by a barrier synchronization. Moreover, the data can be exchanged between threads using a shared memory. The size of a block is limited by the capacity of resources accessible on a single processor core. On currently available GPU’s a single block can contain up to 1024 threads. Each thread block within a grid is uniquely identified by its block ID, and each thread within a block – by its thread ID.
A grid is an array of thread blocks executing the same kernel. The number of blocks in a grid is specified by the amount of data to be processed and the number of available processors. The exchange of data between threads within a grid requires the usage of the global memory. Moreover, while all threads within a single block run simultaneously, different blocks can be executed in any order.
The physical Nvidia GPU architecture consists of several multithreaded streaming multiprocessors (SM in short). Thread blocks are distributed between SM’s in such a way that all threads within a single block are concurrently executed on the same SM (different blocks may, but not necessarily have to, be executed on the same SM). A streaming multiprocessor organizes threads into so-called warps consisting of 32 threads each. The partition is done according to increasing thread ID. Each SM works utilizing SIMT (single-instruction, multiple-threads) architecture, which means that all threads within a single warp execute one common instruction at a time. Any divergence (e.g. caused by a data processed) leads to a serial execution of a single computation path until possible convergence to the same execution path.
A single SM serves multiple warps. It is equipped with a number of warp schedulers and instruction dispatch units (currently 2 or 4 depending on the device used). A scheduler selects a warp, which is ready to be executed and issues it to the physical cores of GPU. If the currently active warp needs to wait for the memory read/write operation it is replaced by another ready warp. While the replaced warp is waiting for the memory access, other warps perform their computations, therefore the SM is busy as often as possible. The waiting period of a single warp is hidden behind the computations of the others (if there are only enough active warps available) and is not seen outside the SM. The simulation of such a behavior is the main part of our model.
The above mentioned heterogeneous environment is equipped with the host memory managed by CPU and the GPU device memory. The latter is significantly more complex than the former. Due to a necessary compromise between the data transfer/access speed and the possible capacity, the GPU device memory consists of various types of data storage, such as global, constant, local and shared memory.
The global memory is the largest and at the same time the slowest type of GPU memory (with hundreds of cycles latency). Together with the constant memory it is the only type of GPU memory, which can be accessed by the host. It is available for reading and writing for all running threads, however the data exchange and result sharing are possible only after a kernel-wide global synchronization.
The content of the global memory is accessed in blocks of size 32, 64 or 128 bytes (depending on the device used). Every time an element within a block is accessed, the whole block have to be transferred. Therefore, the concurrent (among the threads within a single warp) global memory read and write operations are grouped into transactions, the number of which depends on the cache lines required to serve all threads within a warp. However, if different threads in a warp refer to different memory blocks, all such blocks have to be transferred to cache sequentially. Such a situation cause the necessity of repeated global memory accesses.
The shared memory is a fast memory physically placed inside a multiprocessor. It consists of blocks, each of which is available for all threads within a single thread block. Moreover, each such block is divided into several so-called memory banks. The access of different threads to different banks is realized simultaneously, while the access of different threads to the same bank is realized sequentially. Such a situation is called the bank conflict. The shared memory can be used for data exchange between threads within the same thread block after block-wide thread synchronization.
The local memory of a single thread consists of a number of registers, which are the fastest type of memory available (with almost negligible access time). It is used to store local thread variables. Due to large number of threads the capacity of each thread local memory is strongly limited.
To complete the picture we have to mention also the constant and texture memory – dedicated parts of the GPU device memory (usually buffered). Both of them are optimized for access speed within the device, but are available for threads in read-only mode. Neither of them is considered in our model. 3
The set of non-negative integers is denoted by IN. Given a set X, the cardinality (number of elements) of X is denoted by jXj, the powerset (set of all subsets) by 2X – the cardinality of the powerset is 2jXj. Multisets over X are members of INX , i.e., functions from X into IN. For convenience and readability, if the set X is finite, multisets in INX will be represented by vectors of INjXj (assuming a fixed ordering of the set X). The addition and the partial order on INX are understood componentwise, while < means and 6=.
Let us now recall basic definitions and facts concerning inhibitor Petri nets [1,
Definition 1. An inhibitor1 place/transition net (p=t-net) is a quintuple S = (P; T; W; I; M0), where: – P and T are finite disjoint sets, of places and transitions (actions), respectively; 1 Note that in the case of bounded nets the use of inhibitors is not necessary, one can provide an equivalent (with more complex structure) net without inhibitors. – W : P T [ T P ! IN is an arc weight function; – I P T is an inhibition relation; – M0 2 INP is a multiset of places, named the initial marking.
For all a 2 T we use the following denotations: a = fp 2 P j W (p; a) > 0g – the set of entries to a a = fp 2 P j W (a; p) > 0g – the set of exits from a a = fp 2 P j (p; a) 2 Ig – the set of inhibitor places for a.
Petri nets admit a natural graphical representation. Nodes represent places and transitions, arcs with classical arrow heads represent the weight function, while arcs with small circles as arrowheads represent inhibition relation. Places are indicated by circles, and transitions by boxes. Markings are depicted by tokens inside the circles, the capacity of places is unlimited. However, Petri nets used in our model are bounded (which means that there exist a common bound for all the numbers of tokens appearing during the computation in a single place).
The set of all finite strings of transitions is denoted by T , the empty string is denoted by ", the length of w 2 T is denoted by jwj, number of occurrences of a transition a in a string w is denoted by jwja.
Multisets of places are called markings. In the context of p/t-nets, they are typically represented by nonnegative integer vectors of dimension jP j, assuming that P is totally ordered.
A transition a 2 T is enabled at a marking M whenever a M (all its entries are marked) and 8p2 a M (p) = 0 (all inhibitor places are empty). If a is enabled at M , then it can be executed. A marking M is called a dead marking if no transition is enabled at M (which means that 8a2T 9p2P W (p; a) > M (p) _ ((p; a) 2 I ^ M (p) > 0) ). The execution of an enabled transition a is not forced and changes the current marking M to the new marking M 0 = (M a) + a (tokens are removed from entries, then put to exits). We shall denote M a for the predicate “ a is enabled at M ” and M aM 0 for the predicate “ a is enabled at M and M 0 is the resulting marking”.
In this paper however, we use the maximal concurrent semantics and in every marking we execute one of the maximal sets of enabled transitions (i.e. a step, which is maximally concurrent at this marking). Formally, a set of transitions A T is called step and is enabled if Pa2A a 0. The execution of a step A changes the current marMkinagndM8pto2(tShae2Ane wa) mMa(rpk)in=g M 0 = M Pa2A a + Pa2A a . We say that a step is maximally concurrent at marking M if A is enabled at M and 8a2=AA [ fag is not enabled at M .
The notions of enabledness and execution we extend, in a natural way, to strings of steps (computations ): the empty string " is enabled at any marking, a string w = Av is enabled at a marking M whenever M AM 0 and v is enabled at M 0. The predicates M A, M w, M AM 0 and M wM 0 are defined like for single transitions. 4
The Petri net model of GPU calculations is consistent with the workflow of
multithreaded streaming multiprocessors (SMs). The model represents the way
one SM operates. It models each warp assigned to the SM, selection of the
next instruction for the SM, accesses to the global and shared memory, and
arithmetic operations. It does not represent threads hierarchy (blocks, grid),
repeated accesses to the global memory nor bank conflicts. The model allows
simultaneous accesses to the global memory, but the number of warps, which
can use the global memory, at the same time, is limited by the number of SM
warp schedulers (2 or 4). Each element of the model was created basing on [
The size of the considered Petri net depends on the number of warps. The two elements: the place p0 – SM and the transition t0 – waiting are the constant part of the model, other are generated for every warp. The place p0 represents the streaming multiprocessor and its initial marking should correspond to the number of warp schedulers. For the modern graphical cards it should be 2 or 4. This place is connected by a loop with the transition t0. The transition t0 can be executed only when the warp schedulers cannot schedule any instruction, i.e. there is no instruction ready to be executed. The waiting transition is added to the model to gain control over how many steps of the calculations on the SM is idle. The minimization of that number is crucial for optimization of GPU programming. Beside those two elements, the PN model contains 18 places and 17 transitions for every warp required by the analysed algorithm. That part is called a warp part of PN model (WPNM). The detailed description of the most important places and transitions of WPNM is presented in Table 1.
The WPNM together with SM place and waiting transition are depicted in
Figure 2. The place p1 represents the activation of the warp. It is marked when
the warp is active. The place p2 is marked when the warp finishes the execution
of its previous instruction. The warp can be scheduled for the execution only
when the places p2 and p4 are marked. A token in p4 means that the next
instruction is selected and ready. The places from p3 to p10 and the transitions
t2, t3, t4 (the frame I part in Figure 2) are responsible for controlling and selecting
the instructions. The three types of instructions are allowed in the model: an
arithmetic calculation, an access to the shared memory and an access to the
global memory. The initial marking of the place p5 corresponds to the number
xa of arithmetic operations required by the analyzed algorithm. Similarly, the
initial markings of the places p7 and p9 represent respectively the numbers xs
and xg of read/write operations from/to the shared and global memory. If the
place p3 is marked and at least one of the places: p5, p7, p9 is not empty, the
next instruction can be selected. The selection is random. When the instruction
is chosen, according to its type (arithmetic calculation, shared memory access,
global memory access) the marking of the corresponding place is decreased and
a token is added to p4. Moreover, (according to the instruction type) one of
the places: p6, p8 or p10 is marked. Now the selected instruction is ready to be
executed and the warp can be processed by SM. The execution of the instruction
is represented by the three parts of the net marked in frame II, frame III and
frame IV (see Figure 2). They correspond to the type of the selected instruction:
frame II for an arithmetic calculation, frame III – an access to the shared memory
and frame IV – an access to the global memory. The arithmetic calculation is
simply represented by one place and two transitions. When the calculation is
done, tokens are put in places: p2 and p3 (which means that the instruction is
finished and the next one can be selected), and in the place p0 (which means
that SM is ready to execute the next instruction). The same situation is obtained
when the access to the memory (shared or global) is finished, but those parts
of the model contain more places and transitions. Those additional elements are
used to model the memory access latencies. The shared memory latency and the
global memory latency are the parameters of the model and are denoted by l1
for the global memory and by l2 for the shared memory. The transition t9 (t12
respectively) may be executed only after l1 (l2 respectively) executions of t8 (t17
respectively). For testing, their default values were 20 for l1 and 2 for l2, which
is consistent with [
The transition t16 is connected by inhibitor arcs with the places p5, p7 and p9 (the frame I in Figure 2) and can be executed only when those places are empty, i.e. there is no instruction left for execution. The transition is also connected by a regular arc with p2, Moreover, t16 is the only one able to take the token from the place p1 and its execution is equivalent to the termination of a given warp.
As it was mentioned above, the WPNM (places from p1 to p18 and transitions from t1 to t17) is generated for a single warp. In the case of multiple warps a separate WPNM should be generated for each of them. To distinguish between distinct WPNMs one can either increase the numeration of places and transitions accordingly or assign to them two-part labels consisting of the original place/transition number together with the warp id. It should be noticed that each place corresponding to p2 and representing the readiness of the given warp should be connected by an inhibitor arc with the transition t0. Moreover, each transition corresponding to t1 should be connected with the place p0 (SM ). The same goes for transitions corresponding to t5, t14 and t15. Note that the control is returned by t5 not t11 in the case of part responsible for the access to the global memory. Thanks to that, SM can process the next ready warp while the current one is waiting for the memory access.
Our PN model of GPU computation and memory access may by adapted
for any algorithm. Instantiations of the model for different kernels may differ
in the number of warps and the marking of places responsible for instructions
counting (i.e. p5, p7 and p9). For the chosen number of required warps, a new
model containing sufficient number of WPNMs can be generated. The other
possibility is the generation of one big model with the maximal possible number
of WPNMs (i.e. 64 for modern graphical cards [
t12 p16 2 l2 t15 t17 t14 (III) different types of instructions: arithmetic calculation, access to the shared memory and access to the global memory (respectively).
p0 Streaming Multiprocessor (SM) p1 Active warp p2 The previous instruction is finished and
the warp is ready for the next one p3 Check the next instruction p4
Next instruction is ready p5 Arithmetic operations p6 Instruction – arithmetic calculation p7 Access to the shared memory p8 Instruction – shared memory access p9 Access to the global memory p10 Instruction – global memory access that the maximal number of WPNMs is limited to 64 in modern GPUs. The number of arithmetic operations and accesses to the shared and global memory need to be calculated for the considered algorithm. Those numbers should be used as the initial marking of places p5, p7 and p9. If the model is constructed according to the description above, it is ready to be used out of the box. The usage of the model involves the execution of computations according to the maximal concurrent semantics (i.e. concurrent execution of all transitions that are enabled) starting from the initial marking until reaching the dead marking. The latter is obtained only when all places corresponding to p1 (for each WPNM) become empty, which is equivalent to the termination of all warps. The number of steps of the computation is returned by the model and corresponds to the GPU execution time.
The maximal concurrent semantics requires execution of all enabled transitions. However, some of the enabled transitions may be in a conflict (i.e. execution of one transition makes disables another transition). In our implementation of the model, in every step, a permutation of the enabled transitions is randomly generated (using uniform distribution). Transitions are executed according to an order determined by the generated permutation. If two (or more) transitions are in a conflict, the transition appearing earlier in the considered order is executed.
The initial tests of the PN model were performed using the matrix
multiplication kernel from [
In the main experiment we used a binary version of the least significant digit
radix sort algorithm [
In the parallel version we made n runs (one for every bit), each run consisting of three phases. At the beginning of each run, we partition the dataset equally between m nodes. During the first phase, j-th node counts zeros[i; j] – the number of elements containing 0 on i-th bit (consequently, we know ones[i; j] – the number of elements with 1 on i-th bit for this node). In the next phase, we need to compute the positions for the set of elements assigned to each node. Namely, j-th node should place all the elements having 0 on i-th bit between Pk<j zeros[i; k] and Pk j zeros[i; k] 1, while those with 1 on i-th in the The last, third phase, is the rearrangement of the list of integers. Each node traverse assigned part of the data splitting it into two parts (containing only 0 on i-th bit and only 1 on i-th bit) with the use of the positions computed in the second phase and in a stable manner. Since the output space for the nodes is partitioned into disjoint blocks, this phase may be realized using either shared or global memory.
We consider three CUDA implementations of the algorithm described above. In all versions, the array of integers A to be sorted is stored in the global memory. During each kernel execution, one block of threads (with different number of threads) is created. Each thread has its own part of the array A assigned. Its size is the parameter and is denoted by memsize. The product: threadsN umber memsize should be equal to the size of A. At the beginning of each run, every thread copies the assigned part of the array A from the global memory to its local registers, then calculates number of 0 and 1 bits. At the end of the run, the content of the global memory array is rearranged – each thread moves the elements from its part of A.
Each of the three implementations of the radix sort algorithm was tested on a randomly generated array of 65536 integers, with five combinations of threadsN umber and memsize parameters. For the given implementation and the value of parameters, the numbers of arithmetic operations and accesses to the global and shared memory were calculated and used in the PN model (as the initial marking of the places p5, p7 and p9). As an arithmetic operation we count every assignment, addition, subtraction, multiplication, division, relational operation, logical operation and array subscript (arrays in registers). Every access (read or write) to data stored in the global or shared memory is counted as single memory operation.
The numbers of steps of the model calculations were compared to the execution times of kernels. The tests were performed on NVIDIA GeForce GTX 960M graphical card with CUDA Toolkit 8.0. The execution times of kernels are averages of one hundred runs. The results for the PN model were calculated as averages of ten computations of the models.
In the first implementation only the global memory is used. The second phase of the algorithm is performed by thread with id 0. The results of the tests are depicted in Figure 4 – the dotted line.
) s m ( e m i t e g a r e v a l e d o m N P f o s p e t s 60 40 20 0 8 6 4 2
Implementation I Implementation II
Implementation III 64
128 106
Implementation I Implementation II Implementation III
256 number of threads 512 1024 (a) (b) 64 128
256 number of threads 512 1024
In the second implementation the realization of the second phase is organized in a more efficient way. Instead of computing all sums incrementally, we compute all partial sums (for indexes between p 2q and (p + 1) 2q, where p and q are suitable non-negative integers) and use them as input components for other sums.
Having m = 2r nodes we can compute zeros[i; j] and ones[i; j] for all j m in r + 1 cycles with full system load (using all nodes in every cycle). To do it, we compute in c th cycle specific partial sums of lengths between 2c 1 and 2c 1, an example for r = 2 is depicted on Figure 5. In this sample case z[x::y] denotes Px t y zeros[i; t], while o[x::y] – Px t y ones[i; t], each row corresponds to a single element in table zeros or ones, while in subsequent columns the values of partial sums stored in those elements are given. Each arc between x th and y th row denotes the addition of value kept in x th element to the value kept in y th element, the result is stored in y th element.
More specifically, in the first cycle we compute zeros[i; 2k] + zeros[i; 2k + 1] and ones[i; 2k] + ones[i; 2k + 1] placing results in zeros[i; 2k + 1] and ones[i; 2k + 1] respectively. The number of all operations made in this cycle (the number of arcs between first and second column in example) is m=2 + m=2 = m, hence we can utilize all available nodes to do it at once.
In subsequent cycles we compute longer partial sums using already precomputed ones. This way in c th cycle we compute
X zeros[i; 2ck + t] and
X ones[i; 2ck + t]; 0 t u 0 t u where 2c 1 u < 2c 1, while 0 k < 2r c 1, and store the results in zeros[i; 2ck +t] and ones[i; 2ck +t], respectively. Since after the previous cycle all values P0 t u zeros[i; 2ck +t] and P0 t u ones[i; 2ck +t] are stored in memory, where 2c 2 u < 2c 1, while 0 k < 2r c+1, we need to add only two elements for each longer partial sum computed in c th cycle. Note that the number of such operations equals to the size of the range of u multiplied by the size of the range of k and doubled (we need to compute both ones and zeros), i.e. jfu; 2c 1 u < 2cgj jfk; 0
k < 2r cgj 2 = 2c 1 2r c 2 = 2r = m: FPi ntalklyo,niens[it;hte] floarsteacychclek we ma d.dThtheeexceocmuptiuotnedtimsoesfaorf Pkekrnemlszearnods[tih; ek]retsoulatlsl from the PN model are depicted in Figure 4 – the solid line.
In the third version of the implementation the arrays zeros and ones are stored in the shared memory instead of global. The results of those tests are also presented in Figure 4 – the dashed line.
The results of the PN model calculations for the first implementation (Figure 4 (a)) clearly show that this parallelization is not very efficient as compared to the others, especially for larger numbers of threads. This is confirmed by the execution times of kernels (Figure 4 (b)). One can easily observe that in both plots the number of PN model steps and execution times of kernels for this implementation initially decrease with the increasing number of threads, however z[0::0] z[0::1] z[2::2] z[2::3] o[0::0] o[0::1] o[2::2] o[2::3] + + + + z[0::0] z[0::1] z[0::2] z[0::3] o[0::0] o[0::1] o[0::2] o[0::3] + + + + z[0::0] z[0::1] z[0::2] z[0::3] o[0::0]+z[0::3] o[0::1]+z[0::3] o[0::2]+z[0::3] o[0::3]+z[0::3] Fig. 5. The organization of the second phase of the algorithm for m = 4 (as m = 2r; r = 2) nodes in r + 1 = 3 cycles. for more than 256 threads both of them increase. The similar situation can be also noticed for other implementations, but here the growth is more significant. It can also be observed that for the largest amounts of threads the number of PN model steps increases faster than execution times of kernels. In this case the general direction of changes predicted by our model is consistent with the kernels executions (despite of the lack of the exact match of the plots). The model predicted correctly the most efficient choice of the number of threads, which in this case is 256.
The predictions of the PN model for the 2nd and 3rd implementations are more consistent with the execution times of kernels. In both cases differences between implementations are very small. It is probably caused by a relatively small number of shared memory operations in comparison to accesses to the global memory. For the larger number of threads, the increase in execution times of kernels is more significant than in the results from the PN model. The reliable explanation of such a difference is an overhead for communication between threads. It is clear that such overhead will not be observed in the PN model. However, the general characterization of the results is the same both for the model and the kernels. The PN model predicted correctly also the most effective choice of the number of threads for both implementations, which is 512 threads. 6
The purpose of our PN based GPU computations and memory access model is to help in the analysis and optimization parallel algorithms, which are designed to be implemented on CUDA graphical cards. We do not require the source code to be given as an input, however the algorithm description should be detailed enough to estimate the number of arithmetic operations and accesses to the global and shared memory. Note that the other tools (e.g. the roofline model) also require those information. The expected speedup of the computation from a parallelization can be predicted and compared to other algorithms, even without using of any physical GPU device. The model can help to predict which algorithm is the fastest, how much its modifications can affect the speedup of the computation and whether they are significant enough to include them in the source code.
Any inaccurate results demonstrated by our model might be interpreted as a premise that the algorithm should be improved. As an example recall the presented results for the radix sort algorithm. The predictions of the PN model for the first implementation were not satisfying, and the model clearly showed that a different organization of the second phase of the algorithm improved the time of the computation. On the other hand, using the shared memory in this case was not so beneficial. That was confirmed by the GPU kernels execution times.
Another important advantage of the PN model is the possibility to check how different values of parameters and the level of parallelization can affect the final efficiency. As it can be observed in Figure 4, it is not only a theoretical discussion, the problem is substantial and can result in very different execution times of kernels. With the appropriate number of threads, the GPU calculations were even three times faster. It should be also noticed that for all three implementations presented above, different numbers of threads were the most efficient, hence the selection of the one, universal number of threads is not possible. Our model easily allows to check different number of threads (warps) and its predictions seem to be very accurate. Various values of parameters may also result in different numbers of arithmetic operations and accesses to the memory, and it can be also easily introduced and analyzed by the model.
Our model can freely swap instructions of various types. It allows to check
whether different order of the instructions may improve the algorithm efficiency,
for example by allowing to hide thread waiting periods in calculations. If the
results of the PN model are significantly better than the results from GPU
kernels execution, that possibility should be considered. Naturally, the swap of
the instructions is not always possible because of the nature of calculations. One
of the most important improvements of the model would be the introduction
of a partial order over the set of instructions. This can be achieved by defining
dependence of instructions basing on the access to the same variable (see [
In the model, the bounded inhibitor Petri nets are used. In the future we would like also to consider other PN semantics like colored PNs, timed PNs, etc.. Their nature seems to be suitable for our model and they may help to simplify it or introduce new features.
The PN results are expressed as the number of steps performed, while for the GPU computations we use the execution times of kernels in milliseconds. To compare them directly the special coefficient is required to align one result with the other. However, different nature of various algorithms as well as the lack of research on atomic and comparable in terms of time consumption operations makes the issue of finding the universal coefficient a very hard task.
Although designing of the efficient sorting algorithm was not the aim of
this paper, we described the process of improving the parallel version of the
considered radix sort algorithm. Nevertheless, it is worth to note that providing
its further improvements is possible. The radix sort is quite popular, both in the
most significant digit (MSD), normally together with merge sort subroutine [
Note that for different GPU devices the execution time of a fixed kernel may differ, while the number of steps performed by the model stays the same. It would be useful to prepare a set of benchmarks, which for a given device compute the universal scaling coefficient for this device and the model.