Implementing P Systems Parallelism by Means of GPUs

Implementing P Systems Parallelism by Means of GPUs Jose M. Cecilia2,Jos´eM.Garc´ıa2,Gin´es D. Guerrero2, Miguel A. Mart´ınez–del–Amor1, Ignacio P´erez–Hurtado1, and Mario J. P´erez–Jim´enez1 1Research Group on Natural Computing Department of Computer Science and Artificial Intelligence University of Sevilla Avda. Reina Mercedes s/n, 41012 Sevilla, Spain {mdelamor,perezh,marper}@us.es 2Grupo de Arquitectura y Computaci´on Paralela Dpto. Ingenier´ıa y Tecnolog´ıa de Computadores Universidad de Murcia Campus de Espinardo, 30100 Murcia, Spain {chema,jmgarcia,gines.guerrero}@ditec.um.es Abstract. Software development for Membrane Computing is growing up yielding new applications. Nowadays, the efficiency of P systems simulators have become a critical point when working with instances of large size. The newest generation of GPUs (Graphics Processing Units) provide a massively parallel framework to compute general purpose computations. We present GPUs as an alternative to obtain better performance in the simulation of P systems and we illustrate it by giving a solution to the N-Queens problem as an example. 1 Introduction Membrane Computing is an emerging branch within Natural Computing that was introduced by Gh. P˘aun [24]. The main idea is to consider biochemical processes taking place inside living cells from a computational point of view, in a way that gives us a new nondeterministic model of computation by using cellular machines. Up to now, it has not been possible to have implementations neither in vivo nor in vitro of P systems, so the computation and analysis of these devices are performed by simulators. Therefore, P systems simulators are tools that help the researchers to extract results from models. Since P systems was presented, many software applications have been produced [11]. These simulators have to be as much efficient as possible when handling large problem sizes. Thus, the massively parallel nature of P systems computations points out to look for a massively parallel technology where the simulator can run efficiently. Parallel computation on clusters is the traditional environment to speedup parallel applications. Particularly, many simulators of P systems have been designed for clusters of computers [4]. However, this computation is relatively expensive and it is available for organizations that have enough resources to buy and maintain those clusters. Nowadays, there are other cheaper solutions in the computer market that also provides parallel environments. Among these solutions, the newest generation of graphics processor units (GPUs) are massively parallel processors which allow to develop a wide range of parallel applications. We also recall that other parallel computing platforms for P systems simulators are being investigated, such as special hardware circuits [6] and FPGAs [20]. GPUs can support several thousand of concurrent threads providing a massively parallel environment where parallel applications can obtain huge performance [14][17][29]. Current Nvidia’s GPUs, for example, contain up to 240 scalar processing elements per chip [16], they are programmed using C and CUDA [32][21], and they have low cost compared with a cluster of computers. In this paper, we use CUDA as parallel programming environment for P systems simulator in order to speedup the simulation. The input of the simulator is a P system which is defined by using the P-Lingua [5] programming language, and the output is a detailed list of information of every configuration of the computation. The simulation is divided in two main stages: selection stage and execution stage. At this stage of development, the simulator simulates recognizer P systems with active membranes, the selection stage is executed on the GPU and the execution stage is executed on the CPU. The rest of the paper is structured as follows. In Section 2 several definitions and concepts are given for a correct understanding of the paper. Section 3 introduces the Compute Unified Device Architecture (CUDA) and some concepts of programming on GPUs are specified. In Section 4 we explain the design of the simulator. In Section 5 we implement a solution to the N-Queens problem using the simulator and P-Lingua. Finally, in Section 6 we show some results and compare them with the sequential version of the simulator. The paper ends with some conclusions and ideas for future work in Section 7. 2 Preliminaries Polynomial time solutions to NP-complete problems in Membrane Computing are achieved by trading time for space. This is inspired by the capability of cells to produce an exponential number of new membranes in polynomial time. There are many ways a living cell can produce new membranes: mitosis (cell division), autopoiesis (membrane creation), gemmation, etc. Following these inspirations a number of different models of P systems has arisen, and many of them proved to be computational completeness (they are equivalent in power to Turing machines). In this paper we focus on the model of P systems with active membranes.It is one of the most studied models in Membrane Computing and one of the first models presented by Gh. P˘aun [25]. P systems with active membranes is formed by a membrane structure, where a label and a polarization is associated to each membrane. In this model, every elementary membrane is able to divide itself by reproducing its content into a new membrane. Here we provide a short recall of its features (see [25] for details). The model of P system with active membranes is a construct of the form Π=(O, H, μ, ω1,...,ω m,R), where m≥1 is the initial degree of the system; O is the alphabet of objects,His a finite set of labels for membranes; μis a membrane structure (a rooted tree), consisting of mmembranes injectively labelled with elements of H,ω1,...,ω mare strings over O, describing the multisets of objects placed in the mregions of μ;andRis a finite set of rules,whereeach rule is of one of the following forms: (a) [a→v]α hwhere h∈H,α∈{+,−,0}(electrical charges), a∈Oand vis a string over Odescribing a multiset of objects associated with membranes and depending on the label and the charge of the membranes (evolution rules). (b) a[] α h→[b]β hwhere h∈H,α, β ∈{+,−,0},a, b ∈O(send-in communication rules). An object is introduced in the membrane, possibly modified, and the initial charge αis changed to β. (c) [a]α h→[] β hbwhere h∈H,α, β ∈{+,−,0},a, b ∈O(send-out communication rules). An object is sent out of the membrane, possibly modified, and the initial charge αis changed to β. (d) [a]α h→bwhere h∈H,α∈{+,−,0},a, b ∈O(dissolution rules). A membrane with a specific charge is dissolved in reaction with a (possibly modified) object. (e) [a]α h→[b]β h[c]γ hwhere h∈H,α, β, γ ∈{+,−,0},a, b, c ∈O(division rules). A membrane is divided into two membranes. The objects inside the membrane are replicated, except for a, that may be modified in each membrane. Rules are applied according to the following principles: –All the elements which are not involved in any of the operations to be applied remain unchanged. –Rules associated with label hare used for all membranes with this label, no matter whether the membrane is an initial one or whether it was generated by division during the computation. –Rules from (a) to (e) are used as usual in the framework of membrane computing, i.e., in a maximal parallel way. In one step, each object in a membrane can only be used by at most one rule (non-deterministically chosen), but any object which can evolve by a rule must do it (with the restrictions indicated below). –Rules (b) to (e) cannot be applied simultaneously in a membrane in one computation step. –An object ain a membrane labelled with handwithchargeαcan trigger a division, yielding two membranes with label h, one of them having charge β and the other one having charge γ.Notethatallthecontentspresentbefore the division, except for object a, can be the subject of rules in parallel with the division. In this case we consider that in a single step two processes take place: “first” the contents are affected by the rules applied to them, and “after that” the results are replicated into the two new membranes. –If a membrane is dissolved, its content (multiset and interior membranes) becomes part of the immediately external one. The skin is never dissolved neither divided. Note that P systems can be seen as devices with two levels of parallelism: among membranes (every membrane works independently, with the exception of when there are communication across them) and among objects inside a membrane (the rules are applied to the existing multiset of objects in a maximal parallel way). Recognizer P systems were introduced in [26], and constitute the natural framework to study the solvability of decision problems. The data representing an instance of the problem has to be provided to the P system to compute the appropriate answer. This is done by codifying each instance as a multiset placed in an input membrane. The output of the computation, yes or no,issenttothe environment in every halting configuration. Furthermore, the act of simulating something generally entails representing certain key characteristics or behaviours of some physical, or abstract, system. However, an emulation tool duplicates the functions of one system by using a different system, so that the second system behaves like (and appears to be) the first system. With the current technology, we can not emulate the functionality of a cellular machine by using a conventional computer to solve NP-complete problems in polynomial time, but we can simulate these cellular machines, not necessarily in polynomial time, in order to aid researchers. However, depending on the underlying technology where the simulator is executed, the simulations can take too much time. The technology used for this work is called CUDA (Compute Unified Device Architecture). CUDA is a co-designed hardware and software solution to make easier developing general-purpose applications on the Graphics Processor Unit (GPU) [34]. GPUs, that are one of the main components of traditional computers, originally were specialized for math-intensive, highly parallel computation which is the nature of graphics applications. These characteristics of the GPU were very attractive to accelerate scientific applications which have massively parallel computations. However, the problem was the way to program general purpose applications on the GPU. This way involved to deal with GPUs designed for video games, so they have had to tune their applications using programming idioms tied to computer graphics, programing environment tightly constrained, etc [17] [14]. The CUDA extensions developed by Nvidia provides an easier environment to program general-purpose applications onto the GPU, because it is based on ANSI C, supported by several keywords and constructs. ANSI C is the standard published by the American National Standards Institute (ANSI) for the C programming language, which is one of the most used. P systems devices are massively parallel, what fits, in a similar way, into massively parallel nature of the GPUs with thousands of threads running in parallel. These threads are units of execution which execute the same code concurrently on different pieces of data. 3 Graphics Processing Unit Driven by the video games market, programmable GPUs (Graphics Processing Units) have evolved into a highly parallel, multithreaded, manycore processor. They were designed to accelerate graphics applications, which transform threedimensional data (coordinates of triangle vertices) into pixels that are displayed on a screen, using for this task programming interfaces such as OpenGL and DirectX. The massively parallel nature of graphics applications and its arithmetic intensity leads the researches to explore more general non-graphics applications onto the GPU, creating a new programming field called GPGPU (General-Purpose on GPUs). GPUs have become an inexpensive and readily available single-chip massively parallel system. However, GPGPU programmers had to deal with the limitations and difficulties of constrained graphics primitives to compute their non-graphics computations. The emergence of Compute Unified Device Architecture (CUDA) [34] programming model, proposed by Nvidia Corporation in 2007, has helped to develop highly-parallel applications onto the GPU easier than it was before. CUDA allows GPGPU programmers to develop their applications in a more familiar environment by using C/C++ programming language, with some extensions to manipulate special aspects of the GPU. Moreover, Nvidia consolidated this trend launching a line of GPUs optimized for general purpose computations called TESLA [16]. In this work we use a Tesla C1060 graphics processor unit (GPU) from Nvidia as hardware target for its study. This section introduces the Tesla C1060 computing architecture. In addition, it analyses the threading model of Tesla architectures, and also the most important issues in the CUDA programming environment. 3.1 Tesla C1060 Base Microarchitecture The Tesla C1060 [16] is based on a scalable processor array which has 240 streaming-processor (SP) cores organised as 30 streaming multiprocessor (SM). The applications start at the host side (the CPU) which communicates with the device side (the GPU) through a PCI-Express x16 bus (see the top of figure 1). The SM is the processing unit, and it is unified graphics and computing multiprocessor. Every SM contains eight SPs arithmetic cores, one double precision unit, 16-Kbyte read/write shared memory, a set of 16384 registers, and access to the off-chip memory (global/local memory). The access to shared memory is very cheap, however, the access to the off-chip memory has low performance because it is out of the chip, as it is shown on figure 1. In addition, table 1 shows all memories available on the GPU and also the cost to access them. Fig. 1. Tesla C1060 GPU with 240 SPs: Streamming Processors, organised in 30 SMs: Streamming Multiprocessors Table 1. Memory System on the Tesla C1060 Memory Location Size Latency Access Registers On-Chip 16384 32-bits Registers per SM ≃0cycles R/W Shared Memory On-Chip 16 KB per SM ≃registers R/W Constant On-Chip 64 KB ≃registers R Texture On-Chip Up to Global >100 cycles R Local Off-Chip 4GB 400-600 cycles R/W Global Off-Chip 4GB 400-600 cycles R/W 3.2 Parallel Computing with CUDA The GPU is seen as a cooprocessor that executes data-parallel kernel functions. The user creates a program encompassing CPU code (Host code) and GPU code (Kernel code). They are separated and compiled by nvcc (Nvidia’s compiler for CUDA code) as shown in figure 2. Firstly, the host code is responsible for transfering data from the main memory (RAM or host memory) to the GPU memory (device memory), using CUDA instructions, such as cudamemcpy. Moreover, the host code has to state the number of threads executing the kernel function and the organization of them. Threads execute the kernel code, and they are organized into a three-level hierarchy as it is shown in figure 3. At the highest level, each kernel creates a single grid that consists of many thread blocks. Each thread block can contain up to 512 C/C++ with CUDA Extensions NVCC CPU Code PTX Code PTX to Target Compiler G80 PTX Code T10 Fig. 2. Nvcc compilation process threads, which can share data through Shared Memory and can perform barrier synchronization by invoking the –syncthreads primitive [31]. Besides, thread blocks can not perform synchronization. The synchronization across blocks can only be obtained by terminating the kernel. Furthermore, the host code calls the kernel function like a C function by passing parameters if it is needed, and also by specifying the number of threads per block and the number of blocks making up the grid. Each block within the grid has their own identifier [22]. This identifier can be one, two or three Fig. 3. Thread organization in CUDA programming model dimensions depending on how the programmer has declared the grid, accessed via .x,.y,and.z index fields. Each thread within the block have their own identifier which can be one, two or three dimensions as well. Combining thread and block identifiers, the threads can access to different data address, and also select the work that they have to do. The kernel code is specified through the key word global and the syntax is: global kernelName <<< dimGrid, dimBlock >>> (...parameter list...) where dimGrid and dimBlock are three-elements vectors that specify the dimensions of the grid in blocks and the dimensions of the blocks in threads, respectively [21]. 3.3 Threading Model A SM is a hardware device specifically designed with multithreaded capabilities. Each SM manages and executes up to 1024 threads in hardware with zero scheduling overhead. Each thread has its own thread execution state and can execute an independent code path. The SMs execute threads in a Single- Instruction Multiple-Thread (SIMT) fashion [16]. Basically, in the SIMT model all the threads execute the same instruction on different piece of data. The SMs create, manage, schedule and execute threads in groups of 32 threads. This set of 32 threads is called Warp. Each SM can handle up to 32 Warps (1024 threads in total, see table 2). Individual threads of the same Warp must be of the same type and start together at the same program address, but they are free to branch and execute independently. Table 2. Major Hardware and Software Limitations programing on CUDA Configuration Parameters Limitation Threads/SM 1024 Thread Blocks/SM 8 32-bit Registers/SM 16384 Shared Memory/SM 16KB Threads/Block 512 Threads/Warp 32 Warps/SM 32 The execution flow begins with a set of Warps ready to be selected. The instruction unit selects one of them, which is ready for issue and executing instructions. The SM maps all the threads in an active Warp per SP core, and each thread executes independently with its own instructions and register state. Some threads of the active Warp can be inactive due to branching or predication, and it is also another critical point in the optimisation process. The maximum performance is achieved when all the threads in an active Warp takes the same path (the same execution flow). If the threads of a Warp diverge, the Warp serially executes each branch path taken, disabling threads that are not on that path, and when all the paths complete, the threads reconverge to the original execution path. 4 Design of the Simulator for Recognizer P Systems In this section we briefly describe the simulator of recognizer P systems with active membranes, elementary division and polarization. Firstly, we explain the previous work that we have done in order to prepare the development of the parallel simulator on the GPU. Then, we introduce the algorithm design in the CUDA programming language, and finally, we finish with our simulator’s design. 4.1 Design of the Baseline Simulator As previously mentioned, CUDA programming model is based on C/C++ language. Therefore, the first recommended step when developing applications in CUDA is to start from a baseline algorithm written in C++, where some parts can be susceptible to be parallelized on the GPU. In this work, we have based on the simulator for P systems with active membranes developed in PLinguaCore [5]. This sequential (or single-threaded) simulator is programmed in JAVA, so the first step was to translate the code to C++. The simulator is executed into two main stages: selection stage and execution stage.Theselection stage consists of the search for the rules to be executed in each membrane. Once the rules have been selected, the execution stage consists of the execution of these rules. The input data for the selection stage consists of the description of the membranes with their multisets (strings over the working alphabet O, labels associated with the membrane in H,etc.),andthesetofrulesRto be selected. The output data of this stage is the set of selected rules. Only the execution stage changes the information of the configuration. It is the reason because execution stage needs synchronization when accessing to the membrane structure and the multisets. At this point of implementation, we have parallelized the selection stage on the GPU, and the execution stage is still executed on the CPU because of the synchronization problem. We also have developed an adapted sequential simulator for the CPU (called fast sequential simulator), which has the same constraints as the CUDA simulator explained in the next subsections to make a fair comparison among them. This simulator achieves much better performance than the original sequential simulator. 4.2 Algorithm Design in CUDA Whenever we design algorithms in the CUDA programming model, our main effort is dividing the required work into processing pieces, which have to be processed by TB thread blocks of Tthreads each. Using a thread block size of T=256, it is empirically determined to obtain the overall best performance on the Tesla C1060 [28]. Each thread block access to one different set of input data, and assigns a single or small constant number of input elements to each thread. Each thread block can be considered independent to the other, and it is at this level at which internal communication (among threads) is cheap using explicit