Controllers: An abstraction to ease the use of hardware accelerators

3.1 Kernel management

Kernel launching is an operation that inserts an order to execute a kernel in the controller queue, with given real parameters, when the associated device (accelerator or set of CPU cores) is available. This is done using a launching primitive (see Figure 1). The controller internals will ensure that the input data has been transferred to the device memory if needed, and that previous kernels have ended, before doing the real launch. The controller execution policy could reorder the kernels in the queue to maximize the execution efficiency as long as the input/output dependencies across kernels are not violated. The main features and contributions of our model in the kernel management are described below.

3.2 Data management

Accelerators may have their own memory spaces, forcing the transfer of the input data of the kernels, and the obtained results, between the memory of the host platform and the memory of the device accelerator. The manual management of these data movements can be cumbersome and error-prone. Moreover, it is difficult to predict in advance when asynchronous data transfers are possible, or when data should stay in the accelerator device memory, as this depends on the exact sequence of kernels launching. Our model abstracts all these issues from the programmer. The main features/contributions of our model on data management are the following.

4.1 Data structures and Hitmap

Regarding the data structures that the model handles, we have decided to integrate our implementation with Hitmap (Gonzalez–Escribano et al., 2014), an efficient library for hierarchical tiling and mapping of arrays. Hitmap defines the HitTile structure, an abstract entity for n-dimensional arrays and tiles. A HitTile structure is a handler to store array meta-data, along with the pointer to the actual memory space of the data. The circles for the variables A, B, and C in Figure 1 represent the handlers. The HitTile structure should be specialized for each array base type at the beginning of the program.

There are only four functions of Hitmap needed to work with the controllers implementation. First, hit_tileDomain and hit_tileDomainAlloc are used to declare the index domains of a tile array, also allocating the memory for the data in the second case. The function hit_tileFree is used to free the data memory and clean the handler. The function hit_tileElem is used in host or kernels code to access the elements of a tile. It receives a tile name, a number of dimensions, and the index values of the desired element. This function provides a homogeneous interface to manage data structures in both host and accelerator device kernels. The data are stored and accessed in row major order in all cases.

Hitmap library includes many other functionalities. They include management of hierarchical subselections of parts of the tiles, and transparent management of distributed-arrays, with abstract partition and communication functionalities that internally use a message-passing paradigm (exploiting MPI). This will allow the transparent integration of the controllers in distributed multi-node clusters in the future.

Most of the meta-data in the HitTile structure are only needed for advanced functionalities, such as distributed-array partitions and communications in the host code. In our new implementation, we define a new much smaller handler (HitKTile) with the minimum information needed for data accesses to multidimensional arrays through a data pointer. The kernel launching interface will transparently transform the handlers to this new type, substituting the data pointer with its equivalent in the device memory space. The data accesses inside the kernels uses this internal pointer along with a minimum number of arithmetic operations, exposing the expressions to the native compiler to open the possibility of further optimizations. The result obtains as good a performance as the classical array accesses.

In Hitmap, the implementation of different kinds of tiles hides the details of the data access for different internal data layouts to the programmer. The Hitmap library already integrates sparse domains and sparse data structures into the HitTile abstraction. The future transformation of these handler structures, to implement efficient and portable kernels, should follow a similar approach as the one used for dense arrays.

4.2 Controllers and variables management

4.3 Declaration and configuration of kernels

A kernel is declared by using the primitive KERNEL_<type>. Where type may be empty to indicate a generic kernel, usable on any kind of device, or a specific value for a specialized code for a given type of device. This is useful when different optimizations on the kernel code are required for different devices. Currently, the library supports the specific primitives KERNEL_GPU for CUDA code targeting NVIDIA’s GPUs, KERNEL_CPU for host machine code targeting sets of CPU cores, and KERNEL_GPU_WRAPPER for host machine code which includes calls to specialized GPU libraries, like, for example, cuBLAS routines. The controller ensures that both kinds of GPU kernels are executed with exclusive control of the associated device. This allows one to automatically coordinate the launch of sequences of classical kernels and calls to GPU libraries in the same device with our model.

The kernel-definition primitives declare in brackets the number of parameters of the kernel, with a tuple of information for each parameter. The parameter information includes its type, name, and input/output role:

For the case of accelerator devices, with separated memory spaces, this configuration allows the controllers to determine if it is necessary to carry out data transfers with the main memory of the host when kernels are launched. The primitive is followed by a structured block with the kernel code. Thus, it resembles a C function header.

The controller CPU implementation contains the loops to execute the kernel code for each element of the threads index space, that is internally assigned to a specific OpenMP thread. Both CPU and GPU kernels code use a predefined thread structure with three integer elements x,y,z indicating the 3-dimensional indexes of the corresponding fine-grain thread. In order to ease the kernel reuse across different types of device, these indexes are adapted to have the same row-major meaning in both types of kernels. Also, the space of valid indexes that can be defined in the kernel launching is independent of the threadblock sizes in GPUs. Actual threads with identifiers outside the chosen launching index space will be clipped transparently by the controller launching system. For specialized GPU kernels, the kernel code may use CUDA primitives mixed with the Hitmap data accesses that use the new portable index system.

4.4 Kernel characterization

The kernel characterization is a programmer hint to help the system to automatically determine proper kernel launching parameters in terms of special code features and platform architecture information. The CPU thread granularity in our prototype is determined by a simple regular blocking policy that does not require a specific kernel characterization.

For GPU kernels, our current prototype library integrates the model presented in Ortega–Arranz et al. (2014) and Torres et al. (2013). This model allows us to determine configuration parameters (grid, threadblock and L1 cache memory sizes), for NVIDIA’s GPUs. The primitive KERNEL_CHAR, taken from Ortega–Arranz et al. (2015), is used to provide the characterization of the kernels to the controller. The primitive receives the kernel name, the number of dimensions of the thread space (1, 2, or 3), and descriptive values for the characterization model. These values are a qualitative description of characteristics of the kernel code provided by the programmer. They are related to: (a) The coalescing property of the global memory access patterns (full, medium, scatter); (b) The ratio of arithmetic/logic operations per global memory access (high, medium, low); and (c) The ratio of data sharing accesses in a block per global memory access (high, medium, low). For the default case, when the programmer cannot provide a proper characterization for all the parameters, a def keyword can be used instead of one of the given values, and the model provides typical threadblock values for each CUDA architecture, that work well in a general case, maximizing occupancy if possible, etc. For a more detailed description of this qualitative descriptors with tentative quantitative ranges, see the works Ortega–Arranz et al. (2014), Torres et al. (2013) and Ortega–Arranz et al. (2015).

Our implementation extends this characterization with the possibility to specify a fixed size for the threadblock. This is useful for kernels that rely on specific block sizes or geometries to manipulate shared memory (for example the matrix multiplication code in the CUDA Toolkit Samples).

4.5 Kernel launching

The function CtrlLaunch is used to launch a kernel, with given real parameters, to the computational device associated to a controller. The launched kernel will be enqueued, and eventually executed with the corresponding configuration derived from the information provided by the characterization primitives. Currently, the prototype only supports a First-Come-First-Service policy for kernels execution. The launching function has the following parameters:

The dimensions of the threads index space are specified using a CtrlThreads structure that stores up to three integer values representing the cardinality on each dimension. It is equivalent to the dim3 type in CUDA, but it is used for both, GPU or CPU kernels independently. The variable parameters should be tile variables associated to the controller, internal or bounded, or any value of the proper type for the INVAL parameters.

Internally, the execution of a GPU kernel implies: (1) the creation of the small handlers using the original tile handlers information; (2) The use of the characterization model to select the grid and threadblock geometries; and (3) the execution of a wrapper function. This wrapper reorders the thread indexes to be used as in the CPU kernels code to access data elements in row-major order, and ensures that threads outside of the required index space return immediately before executing the kernel code (see an excerpt of this code on the left of Figure 2). For CPUs, the kernel execution internally implies the partition of the index space in coarse blocks. Figure 2 (right) shows an example of a simplified code that performs this partition. The kernel launching (line 18) calls to an inline function that is generated when this kernel implementation for CPU-cores is declared. A final global synchronization is needed to preserve the launching semantic before proceeding to launch another kernel.

OPEN IN VIEWER

Figure 2:

Excerpts of the controller library code generated for kernel deployment/launching on a CUDA capable GPU device (left), and a group of CPU-cores (right). In the CUDA version the block geometry is selected using the kernel characterization provided by the programmer. We also show the part of the wrapper function launched, that reverses the thread indexes before calling the actual kernel code. For the case of the CPU-cores, each OpenMP thread (assigned to a core) executes this code. In this example, the 2-dimensional index space is divided into blocks by rows without balancing the remaining. The code executes the loops that call the kernel code for each index-element mapped to the thread. It simulates the many-thread approach used for GPU programming using coarser-grain OpenMP tasks.

4.6 Programming example

In this section, we present a practical application of the library concepts described in the previous section by using an example. Figure 3 shows the implementation of a Stencil 2-D application. On the left of the figure, at the beginning of the program, we declare the specialized Hitmap array types to be used in the program (see line 3 in Figure 3 left).

OPEN IN VIEWER

Figure 3:

Example of the kernel characterization and definition (left) and main code (right), for a stencil program implementing an iterative Jacobi PDE solver for the Poisson’s heat diffusion equation. The kernels are usable for both CPU and GPU controllers.

Kernel characterization: Two kernels are characterized at lines 6 and 7 of Figure 3 (left). The first kernel characterization uses the default configuration, while the second one looks for a particular configuration for the updateCell kernel. This last characterization is based on the code properties. According to the classifications presented in Ortega–Arranz et al. (2014) and Torres et al. (2013), the global accesses of the program suggest an almost fully coalesced memory access pattern. The ratio of arithmetic operations per data element is low because less arithmetic operations are carried out compared to read operations. The ratio of data shared between threads of the same block is high, because most neighbour values are reused across their computations.

Kernel definition: The kernel definitions start in lines 10 and 18. The kernel definitions include the name of the kernel, the number of parameters, a specification for each parameter of its input/output mode and type, and the parameter name. Both are generic kernel codes valid for any kind of device as the KERNEL primitive points out. The code in the kernel bodies is executed in parallel on the device by as many threads as we choose in the main program. Another kernel definition example can be seen in Figure 4 (right). In this case, we see a specific kernel code for GPUs that calls a function of a specialized library for NVIDIA’s GPUs.

OPEN IN VIEWER

Figure 4:

Characterization of the generic or GPU specialized kernels for the case studies (left). Example of kernel wrapper to execute a specialized GPU library function (right).

Data management: All the accesses to the variables in the generic kernel body (either reads or writes) must be done through the Hitmap functions. For example, line 14 of Figure 3 (left) shows the copy of an element of the matrix src on the matrix dst. It uses the hit_tileElem function to access the data element corresponding to the indexes assigned to the logical thread that executes the kernel. We can see another example on line 24 of Figure 3 (left). The dst matrix element corresponding to the thread indexes is updated using the values of the neighbour cells in the matrix src.

We show on the right of Figure 3 the code of the main function that coordinates the application execution.

Data structures: Data structures are created and initialized in the first part of the main program (see lines 1 to 9 in Figure 3 right). The M matrix is allocated on the host, where it is initialized. Thus, we use the hit_tileDomainAlloc function to create the data structure. However, the Mcopy matrix is an internal variable, and does not need to have allocation in the host memory space. Only a virtual representation is created, without assigning actual memory, by using the hit_tileDomain function.

Controller entity: Lines 13 and 14 of Figure 3 (right) show examples of the controller creation. One associated to the third GPU of the system, and another one associated to the subset of CPU cores with indexes in the range 4 to 7. Once the controller is created, host data structures are bound to the target devices, and internal data structures are created in the target devices through the controller (see lines 17 and 18 of Figure 3 (right)).

Kernel launching: With the data structures in the target device, we can launch the kernels. The kernel would be executed by as many threads as a CtrlThread object specifies. We see in lines 21 and 22 how two CtrlThread objects are created. The first one includes the domain of the whole data structure, and the second one includes the whole data structure without the borders. After that, in lines 24 and 26, the two kernels are launched using their different thread-index domains, by using the different CtrlThread objects. After the desired number of iterations of the kernels, the control of the result data structure is transferred again to the host by unbinding it, and the controller is destroyed (see lines 31 to 33).

To summarize, we observe that the final program presents an organized sequence of programming phases, that leads to a clear structure. Easy generic guidelines to program with this library can be deduced from this simple example.

5.1 Case studies

As case studies we have selected the following programs. All of them work with floating point numbers. From the CUDA Toolkit Samples, we have selected the Black–Scholes program, and two Matrix Multiplication examples, one using GPU shared-memory, and another directly calling the optimized CUBlas routine. We also use the stencil program that is shown in Figure 3. Finally, we have selected a code to independently apply a trivial recurrence equation to elements in a matrix. A quick description and the motivation to choose these examples follows. The characterizations of the kernels using the proposed model, are shown in the left of Figure 4. The parameters where chosen following the guidance presented in Ortega–Arranz et al. (2014) and Torres et al. (2013).

5.2 Development effort and code complexity

The first part of our experimental study evaluates how the use of our proposed model affects the development effort when compared with using the native programming models for the two types of devices considered in the current version of the prototype: CUDA for GPUs, and OpenMP for multi-core CPUs.

We measure three classical development effort and code complexity metrics: COCOMO lines of code, number of tokens, and McCabe’s cyclomatic complexity. The first two ones measure the volume of code that the programmer should develop. The third one measures the rational effort needed to program it in terms of code divergences and potential casuistry that should be considered to develop, test, and debug. The metrics are applied to the part of the code that includes the kernel, the functions invoked by them, and the host coordination code. We ignore input data initialization, error or results checking, performance instrumentation, and writing messages to the standard output. Thus, the considered host code includes the declaration, creation and initialization of controller, data containers, and structures for parallelism coordination, the memory transfer between host and devices, and the kernel launching operations.

Table 1 (left) shows the measurements for these metrics in the baseline versions, using OpenMP for the CPU variant, and CUDA for the GPU variant, and the versions using our controllers library interface. The results indicate that programming with controllers for GPUs generates a significantly lower volume of code, and a reduced cyclomatic complexity, indicating a clearly lower development effort than using CUDA. A closer look at the codes indicates that most of the reduction is found in the host part of the codes, as expected. On the other hand, for this kind of data parallel codes, the comparisons with OpenMP codes show that using the controllers interface reduces the cyclomatic complexity slightly, but increases the code volume.

OPEN IN VIEWER

Table 1:

Measurements of the development effort metrics for the codes of the case studies (left). Comparison in terms of the percentage of words that are common and can be reused, should be deleted, or should be changed, when porting codes between GPU and CPU versions using the native models, or the controllers model (right).

Case study	Version	Lines of code	#Tokens	Cyclomatic complexity
Recurrence equation	CUDA	44	404	5
	Ctrl.GPU	33	315	3
	OpenMP	27	243	4
	Ctrl.CPU	36	339	3
Black–Scholes	CUDA	148	903	8
	Ctrl.GPU	106	704	7
	OpenMP	81	539	7
	Ctrl.CPU	90	707	6
Jacobi solver	CUDA	43	445	8
	Ctrl.GPU	40	371	4
	OpenMP	33	310	6
	Ctrl.CPU	47	456	5
Matrix mult.	CUDA	71	614	5
	Ctrl.GPU	47	429	4
	OpenMP	22	235	4
	Ctrl.CPU	40	389	4

Case study	CUDA → OpenMP	Ctrl.GPUs → Ctrl.CPUs
Recurrence equation	Common 14%	Common 95%
	Delete 29%	Delete 4%
	Change 57%	Change 1%
Black–Scholes	Common 56%	Common 72%
	Delete 32%	Delete 17%
	Change 11%	Change 11%
Jacobi solver	Common 13%	Common 91%
	Delete 42%	Delete 0%
	Change 44%	Change 9%
Matrix multiplication	Common 8%	Common 49%
	Delete 4%	Delete 3%
	Change 88%	Change 48%

The main advantage of considering the CPU cores as an accelerator device in the controller model is found in the portability of code between GPUs and CPUs. Table 1 (right) shows the percentage of words that can be reused, should be deleted, or should be changed to port from a CUDA program to the equivalent OpenMP version; and the same measurements when porting from a controller version for GPUs to the equivalent controller version for CPUs. The results clearly indicate that the portability of the controller versions across different device types is really high in the first three cases, and still significant for the matrix multiplication program, that includes different specialized kernels. In this case the coordination code using controllers is still the same, deriving in a significantly lesser number of words to change than when porting between CUDA and OpenMP.

5.3 Performance study

The second part of our experimental study measures the impact of using our controllers prototype in terms of performance. We have run the GPU implementations in an NVIDIA’s GeForce Titan Black X, with CUDA capability 5.2, installed on a host machine named Hydra, with two CPUs Intel Xeon E5-2609 v3 @1.90GHz, and 64Gb DDR3 main memory. To test the CPU implementation we have used a shared-memory machine with a higher number of cores. It is named Heracles, and it is a Dell PowerEdge R815 server, with 2 AMD Opteron 6376 processors at 2.3 GHz, with 16 cores per processor, and a total of 32 cores. The operating system in both machines is a Linux Centos 7 OS. The programs have been compiled using the CUDA Toolkit 7.5, and GCC 4.8.3. We have used the flags -arch=compute_52, -O3, and the -mtune choice appropriate for each host machine.

We present comparisons of the execution times skipping input data initialization, results checking, and control messages writing. The instrumented code includes controllers creation and variables binding/unbinding operations (which may imply data transfers). In the same way, CUDA data copies between the host and the device memories are also included. In the case of GPU programs, we also measure the execution time accumulated by the kernels launching alone, without memory transfer operations. Due to instabilities in the execution times of the host codes and data transfers, the results always show the mean of the execution time obtained on five repetitions of each test.

For the CPU versions we have selected input sizes with enough load to achieve scalability in our test machine with up to 32 threads. The exact input size parameters chosen are presented in Table 2. The total execution times of the baseline programs and the versions based on controllers are presented in Figure 5, for different number of active threads/cores. The results indicate that our implementation of the controller abstraction for CPUs does not implies significant overheads for the scalability ranges tested.

OPEN IN VIEWER

Table 2:

Input data sizes and number of iterations selected for each case study in the performance experimental study.

Case study		Recurrence	Black–Scholes	Jacobi	Matrix mult.
Number of iterations		500	512	200	-
CPU data size		15,000 × 15,000	1,000,000	10,000 × 10,000	3000 × 3000
GPU data size	S	10,000 × 10,000	1,000,000	1000 × 1000	3008 × 3008
	M	15,000 × 15,000	10,000,000	5000 × 5000	5024 × 5024
	L	20,000 × 20,000	50,000,000	10,000 × 10,000	10,016 × 10,016
	XL	25,000 × 25,000	100,000,000	20,000 × 20,000	15,008 × 15,008

OPEN IN VIEWER

Figure 5:

Execution times (seconds) in Heracles machine of the baseline (Base) and the controller versions for CPU device (Ctrl.CPU) with a variable number of cores.

For a fair comparison, both GPU approaches (CUDA and controller based) use the same values for the configuration of launching parameters. In the case of the Black–Scholes program, the threadblock size value predicted by the characterization model is the same as that found in the original CUDA code. In the matrix multiplication case, it is fixed in both cases to 32 × 32. For the other two case studies, the values predicted by the characterization model are injected manually in the baseline versions. These values have shown to produce at least 96% of the best performance obtained with other threadblock geometries. We compare the performance obtained when launching the baseline and the controllers versions with different input data sizes, selected to produce from very low to significant execution times (from tenths of a second, to more than four seconds). The exact input size parameters are presented in Table 2, with a class name (S, M, L, XL) for easier reference in the following discussion.

The results for the GPU programs are presented in Table 3. We skip to the present times for the matrix multiplication using the cuBLAS library. The memory transfers and launching operations are enclosed into the library routine, that is called by both the CUDA and the wrapper virtual kernel in the controller version. Thus, the performance is the same except for stochastic behaviours that affect both codes equally.

OPEN IN VIEWER

Table 3:

Execution time (seconds) for the case studies versions using CUDA, or controllers for GPUs, with different input sizes. Host time measures the data communications, coordination code, and launching operations. Kernels time measures the actual kernel execution times. Total times are the addition of both times for easy comparison.

		Recurrence		Black–Scholes		Jacobi		Matrix mult.
Size	Measure	CUDA	Ctrl.GPU	CUDA	Ctrl.GPU	CUDA	Ctrl.GPU	CUDA	Ctrl.GPU
S	Host	0.2580	0.2559	0.0084	0.0047	0.0033	0.0023	0.0357	0.0274
	Kernels	0.5592	0.5723	0.0409	0.0454	0.0151	0.0167	0.0663	0.0666
	Total	0.8172	0.8283	0.0493	0.0500	0.0184	0.0191	0.1020	0.0941
M	Host	0.5720	0.5721	0.0617	0.0654	0.0672	0.0676	0.0495	0.0416
	Kernels	1.2753	1.2752	0.3862	0.3951	0.3203	0.3285	0.2878	0.2983
	Total	1.8473	1.8469	0.4479	0.4604	0.3875	0.3960	0.3373	0.3399
L	Host	1.0249	1.0236	0.2955	0.2899	0.1167	0.1168	0.4977	0.4944
	Kernels	2.2502	2.2502	1.9166	1.9221	1.2277	1.2306	2.1012	2.0925
	Total	3.2750	3.2739	2.2121	2.2120	1.3444	1.3474	2.5989	2.5869
XL	Host	1.5909	1.5864	0.5954	0.5695	1.1901	1.1915	0.4191	0.5748
	Kernels	3.4668	3.4669	3.8311	3.8652	4.8493	4.8574	7.0198	7.0203
	Total	5.0577	5.0529	4.4265	4.4346	6.0394	6.0490	7.4389	7.5441

The execution times of the kernels show a little overhead when using the controller interface, due to an extra internal stage for thread index conversion and threads space clipping. Nevertheless, the overhead is always less than 1% of the kernels mean execution time. The recurrence equation kernel, specifically designed to test the efficiency of the data accesses through the small kernel tile handlers, shows the minimum overhead of the four cases.

As expected, the execution times of the host codes, in both CUDA and controller versions, are smaller than the corresponding kernels execution, specially for large input sizes (L and XL), that present significant kernel execution times. These times spent on host code coordination, and memory transfers, are less stable than the kernels execution times. Sometimes, the controller versions show even better results than the equivalent host code in the CUDA versions. The same memory transfers in the CUDA codes are internally executed in the controller versions. Memory transfers have the less predictable times, with a higher variance. The rest of the controller operations are in general light loaded. The critical one is the kernel launching, that involves several checks and handler manipulations. Nevertheless, the overhead of the launching operation has been measured to be less than 1×10⁻⁵ seconds for the GPU implementation, in the machine used for this study. Thus, the accumulated overhead of many launching operations in the Black–Scholes, or the Jacobi programs, is not significant even for the smaller input sizes tested (class S). This overhead could only be clearly noticeable when executing long sequences of kernels with really low computational load, which are not in general appropriate for parallelism exploitation.