The technique that allows to extend GPU capabilities to deal with data volumes that outfit internal GPU's memory capacity is developed. The technique involves loop tiling and data serialization and could be applied to utilize clusters consisting of several GPUs. Applicability criterion is specified. Transforming scheme designed and semiautomatic proof-of-concept software tool are implemented. Conducted an experiment to demonstrate the feasibility of the proposed approach.
Loop parallelization is a long-standing problem of computational programming. Loops give a fair parallelization opportunity for numerous scientific modeling problems that involve numerical methods. Along with spreading of GPGPU technology [1] that allows employment of graphics accelerators for solving computational tasks new challenges arises. As far as GPU is not a standalone device and is managed by a host operating unit, it should be considered within the context of heterogeneous computational platforms. Composing the programs for such the platforms demands knowledge in architecture and specific programming tools. Generally, the concurrent software development passes through the stage of successive implementation that becomes a starting point for further platform-dependent and hardware environment specific implementations.
The existing automatic code parallelizing tools [2, 3, 4] don't account the limited
amount of GPU's on-board memory space while real-life problems demand in huge
amounts of data to be processed. To embrace those cases of massive computational
tasks that involve large amounts of data we propose a technique that provides an ability
to rip the loop and to split the data and calculation operations.
(
for 0 s0 S0 : for sN L(I N , SN ) iN minsN 1 L(I N , SN ), SN : ...
for s0 L(I0, S0) i0 mins0 1 L(I0, S0), S0 :
F i , D.
where Ik , 0 k N , for 0 i I : S(i) is a notion for a sequence of
S(0), S(
for sn L(In, Sn ) in min sn 1 L(In, Sn ), Sn , L(a, b) a 1 0, a mod b ,
b The inner loop is similar to the initial loop but of a diminished scale. Leaving the inner N+1 loops, let's group the first N+1:
N for 0 e Si :
i0 i g(e); Here g() is a mapping that restores the vector of counters i and is constructed the following way: g0(e) e modS0,
gk (e) e gN (e) e
N 0 e Sk .
k 0
N 1 S j , j0
N jk 1
j1 k 1 g j (e) Sl S j , 0 k N ,
l0 j0
Loop (
Let's denote the inner N+1 loops of the cycle (
N for 0 e Si :
i0
serialize(e, inputData, dataPull);
transfer2device(inputData);
kernel(e, inputData, outputData);
transfer2host(outputData);
deserialize(e, outputData, dataPull);
(
Let's consider the node that consists of one multicore CPU and one GPU. Modern GPUs support direct memory access technology thus allowing to proceed data transfer and kernel execution asynchronously. To optimize data exchange process dual buffering is involved. Four buffers at both host and device sides are involved – two for the input and two for the output data exchange. On the host side, calculations proceed in two threads that execute kernel, serialize and deserialize procedures simultaneously. One of them serialize input data and fills the input data buffer, then transmits the buffer to the GPU and launches the kernel, and the second receives output data buffer from GPU and deserialize it. Besides the calculations, GPU carries bidirectional data transfers through the asynchronous data transfer mechanism. Calculations are performed in three stages – initial, cyclic and finalizing.
At the starting point, data buffers are empty, and GPU awaits the data transfer. It doesn't matter what of the threads will carry the initial step as all of the operations are run successively and asynchronous transfer mode is not involved. At the initial step, CPU serializes input data buffers of the first two iterations and transfers the buffer containing the first iteration data to the accelerator.
After the initial step, the cyclic stage starts. The execution flow is shown at the diagram at Fig.1. On the diagram, the iteration's number of which the data is stored in the buffer is given after the buffer's name. One step of the cyclic lap divides into odd and even parts. Both odd and even parts of the first step skip deserialization as the host output buffers are empty yet. At the odd part of the first step, the accelerator-to-host transfer is omitted too. Meanwhile, an accelerator performs calculations over the current buffer, host threads fetch data buffer from the previous step, deserialize penultimate step buffer, send input data buffer for the next step and serialize buffer for the after the next step. In one step two kernel launches are executed. After each of the parts, odd or even, is finished the processes synchronize. Two final steps depend on the actual kernel launches number. If the number of kernel launches is odd, the final step of the cyclic part excludes an even part and does not involve serialization and host-to-accelerator transfer, and the even part of the penultimate step skips serialization. Otherwise, if the number of launches is even the last cyclic loop step is full, but the even part of the final step omits serialization.
The finalizing step deserializes output data buffer transferred at the last cyclic loop step and then fetches and deserializes the final output data buffer consequently finishing the calculations. 4
In this section, we illustrate the application of the proposed approach to matrix multiplication and N-body problems. The time measurements were collected on the hardware system composed of Intel Core i5-3570 CPU (4 cores 3.8Hz) with 16Gb of host memory and NVIDIA Tesla M2050 GPU (3Gb global memory, 384 bits memory bandwidth, connected through PCIe2.0 x8) running Ubuntu 16.04 host OS.
A semi-automatic source-to-source code transformation tool based on the TermWare
rewriting system [7] aiming to assist in constructing a new concurrent program was
implemented. It takes the initial loop marked with pragmas, applies the transformations
(
The algorithm of the initial sequential matrix multiplication program involved threedimensional loop nest. It was transformed using the proposed technique and C-toCUDA compilers PPCG and Par4All. Both of the programs generated by PPCG and P4A showed comparable results. After applying the slicing technique the initial matrices were split into submatrices. The internal loop subdivision parameter S0 was set to 1, the roles of parameters S2 and S1 is adjusting submatrices width. The schema with double data exchange buffering and two CPU threads was used. Even not involving GPU, adjusting the slicing number allowed to reach about 12 times acceleration over the initial loop due to CPU caching. For the GPU implementation, the parameterized PPCG generated kernel was used; the source codes of the constructed matrix multiplication program are available at GitHub [8]. The chart at Fig. 3 shows the constructed program's timings and the timings of the program obtained with PPCG relatively to the matrix dataset size. The relative acceleration of about 430 times in comparison to the sequential program executed on the CPU for the datasets of square matrices of the size of 5000x5000 single precision floating point numbers was reached. It could be seen from the chart that the PPCG generated program has achieved maximum data set size less than 300 Mb that is 10% of GPU's global memory available. This is due to the fact that PPCG is limited to static memory usage, thus blocking to link programs with too large static arrays. It is worth to mention that involving two CPU threads is excessive as the part of serialization and deserialization is negligible compared to GPU kernel computation time, that could be seen from the GPU execution profile given at Fig.2. Thus involving just one thread instead won't decrease performance substantially, however, two concurrent threads are required to avoid gaps in kernel launches and to gain maximum performance from the GPU. Another application investigated was a predictor routine from N-body problem with predictor-corrector time-iterative [9] algorithm. The model of the system consists of a set of particles that interact pairwise PPCG applying caused a slowdown effect and led in about 500 times decrease in performance in comparison to the sequential CPU implementation. As for the constructed program with a self-implemented kernel not involving shared memory usage, the relative CPU to GPU acceleration at the selected data size range reached 13 times. The plot at Fig.4 shows the dependency of successive CPU and transformed GPU programs execution time on the data size that is scaled by the alteration the number of particles N. The timings are measured for the one time-step of the prediction routine. The memory limit was not reached as it would take approximately 30 years to process one time-step of the algorithm for a fully loaded GPU that was used in the experiment, however, the applicability of the approach is confirmed.
Thus the difference in the constructed multiplication and N-body programs consists in the kernel, serializer and deserializer implementation while the control flow structure remains identical.
Fig. 3. The dependency of the execution time on the size of the input data for the concurrent PPCG-generated and constructed matrix multiplication programs.
This paper proposes an approach for semi-automated parallelization of nested loops for graphics processors. The approach is illustrated by the development of CUDA programs for solving matrix multiplication and N-body problems. As a result, a common unified scheme was used to parallelize both multiplication and N-body problems. An assistant semi-automatic code transformation tool was implemented. Further work relates to developing unified methods and tools in loop parallelization.