Population Dynamics P Systems on CUDA

Population Dynamics P Systems on CUDA Miguel A. Mart´ınez-del-Amor1, Ignacio P´erez-Hurtado1, Adolfo Gastalver-Rubio1, Anne C. Elster2,andMarioJ.P´erez-Jim´enez1 1Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, University of Seville, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain {mdelamor,perezh,marper}@us.es,[email protected] 2HPC-Lab, Department of Computer and Information Science, Norwegian University of Science and Technology, Sem Sælands vei 9, NO-7491, Trondheim, Norway [email protected] Abstract. Population Dynamics P systems (PDP systems, in short) provide a new formal bio-inspired modeling framework, which has been successfully used by ecologists. These models are validated using software tools against actual measurements. The goal is to use P systems simulations to adopt a priori management strategies for real ecosystems. Software for PDP systems is still in an early stage. The simulation of PDP systems is both computationally and data intensive for large models. Therefore, the development of efficient simulators is needed for this field. In this paper, we introduce a novel simulator for PDP systems accelerated by the use of the computational power of GPUs. We discuss the implementation of each part of the simulator, and show how to achieve up to a 7x speedup on a NVIDA Tesla C1060 compared to an optimized multicore version on a Intel 4-core i5 Xeon for large systems. Other results and testing methodologies are also included. Keywords: Ecological Modeling, Population Dynamics, Membrane Computing, Parallel Simulation, GPU Computing, CUDA. 1 Introduction Membrane Computing [19] is part of the broader field of natural or bio-inspired computing. The related computational models are called Psystems.Theyare hierarchically distributed models inspired by how membranes compartmentalize living cells into ”protected reactors”. The model consists of simultaneous applications of rules (abstraction of chemical reactions) over multisets of objects (abstraction of chemical compounds) [20]. Membrane Computing covers both the study of the theoretical basis for the models as well as the applications of the model to various fields including computational Systems Biology [7,22], and Ecosystem Dynamics [3,4,8]. Population Dynamics P Systems, or PDP systems, is a P system based framework for modeling population dynamics [3,17] . It enables simultaneous evolution of a high number of species, as well as the management of a large number of auxiliary objects. It also facilitates model development that can be easily interpreted by simulation software. So far, several algorithms have been developed in order to capture the semantics defined by the modeling framework. A comparison on the performance of these algorithms can be found in [9]. These algorithms select rules according to their associated probabilities, while keeping the maximal parallelism semantics of P systems. In [17], a new simulation algorithm is presented, called DCBA (Direct distribution based on Consistent Blocks Algorithm). It overcomes a common problem on the previous algorithms, regarding a distorted selection of rules. Furthermore, the DCBA was initially implemented and validated using the P- Lingua software framework [11,24] which resulted in a simulation Java library (pLinguaCore). However, the simulations were slow (taking hours to run a single simulation), since the pLinguaCore library is not focused on efficiency. A more efficient implementation based on C++ and OpenMP was presented in [16], taking advantage of modern multicore architectures. These preliminary results indicate that the simulation of PDP systems are memory bound. GPU computing [15] has been successfully used to implement other P systems simulators [2,5,6]. By using CUDA [13,23], the simulators are accelerated taking advantage of the many-core GPU architecture. This paper, introduces our new CUDA-optimized PDP systems simulator. We describe how the data structures have been optimized, and how also the code is adapted and restructured in different parts. A performance analysis is also provided, comparing the results with further optimized sequential and multicore versions written in C++/OpenMP. This paper is structured as follows: Section 2 gives an overview of the P systems based framework that our simulator implements. Section 3 explains the simulation algorithm. Section 4 contains some concepts about the CUDA programming model. Section 5 explains the details of the implementation using CUDA. Section 6 shows some performance results by using a random generator of PDP systems. Finally, conclusions and future work are discussed in Section 7. 2 The P Systems Based Framework Definition 1. A PDP system of degree (q,m)&timeT(q,m,T ≥1)isatuple Π=(G, Γ, Σ, T, RE,μ,R,{fr,j :r∈R, 1≤j≤m},{Mij :1≤i≤q, 1≤j≤m}) where: –G=(V,S)is a directed graph. Let V={e1,...,e m}whose elements are called environments; –Γis the working alphabet and ΣΓis an alphabet representing the objects that can be present in the environments; –Tis a natural number that represents the simulation time of the system; –REis a finite set of communication rules between environments of the form r≡(x)ej pr −−−→ (y1)ej1···(yh)ejh,wherex, y1,...,y h∈Σ,(ej,e jl)∈S (1 ≤l≤h)and pris a computable function from {1,...,T}to [0,1] depending on x, j, j1,...,j h. If for any rule, pris the constant function 1, then we can omit it. These functions verify the following: for each ej∈V and x∈Σ, the sum of functions associated with the rules whose left-hand side is (x)ejis the constant function 1. –μis a membrane structure consisting of qmembranes injectively labelled by 1,...,q. The skin membrane is labelled by 1. We also associate electrical charges from the set {0,+,−} with membranes. –Ris a finite set of evolution rules of the form r≡u[v]α i→u[v]α i,where u, v, u,v∈Γ∗,i(1 ≤i≤q),u+v=λand α, α∈{0,+,−}. •If (x)ejis the left-hand side of a rule from RE, then none of the rules of Rhas a left-hand side of the form u[v]α 1, having x∈u. –For each r∈Rand for each j(1 ≤j≤m),fr,j :{1,...,T}−→[0,1] is a computable function verifying the following: for each u, v ∈Γ∗,i(1 ≤i≤q), α, α∈{0,+,−} and j(1 ≤j≤m)the sum of functions associated with jand the rules whose left-hand side is u[v]α iand whose right-hand side has polarization α, is the constant function 1. –For each j(1 ≤j≤m),M1j,...,Mqj are strings over Γ, describing the multisets of objects initially placed in the qregions of μ,within the environment ej. It is usual to manage multisets through strings. A finite multiset m={af(a1) 1,...,a f(ak) k}can be represented by the string af(a1) 1...a f(ak) kover the alphabet {a1,...,a k}. Nevertheless, all permutations of this string precisely identify the same multiset m. Throughout this paper, we speak about “the finite multiset m”wheremis a string, and meaning “the finite multiset represented by the string m”. That is, a system defined as above can be viewed as a set of menvironments e1,...,e minterlinked by the edges from the directed graph G. Each environment ejcontains a P system, Πj=(Γ, μ, R, M1j,...,Mq,j), of degree q,where every rule r∈Rhas a computable function fr,j associated with it. For each environment ej,wedenotebyRΠjthe set of rules with probabilities obtained by coupling each r∈Rwith the corresponding function fr,j . Aconfiguration of the system at any instant tis a tuple of multisets of objects present in the menvironments and at each of the regions of each Πj, together with the polarizations of the membranes in each P system. We assume that all environments are initially empty and that all membranes initially have a neutral polarization. We assume a global clock exists, sychromnizing all membranes and the application of all the rules (from REand from RΠjin all environments). The P system can pass from one configuration to another by using the rules from m j=1 RΠj∪REas follows: at each transition step, the rules to be applied are selected according to the probabilities assigned to them, and all applicable rules are simultaneously applied in a maximal way. An evolution rule r∈R,oftheformu[v]α i→u[v]α i, is applicable to each membrane labelled by i, whose electrical charge is α, and it contains the multiset v, and its father contains the multiset u. When such rule is applied, the objects of the multisets uand vare removed from the father of membrane iand membrane i, respectively. Simultaneously, the objects of the multiset uare added to the father of membrane i, and objects of multisets vare introduced in membrane i. The application also replaces the charge of membrane ito α. Aruler∈RE,oftheform(x)ej pr −−−→ (y1)ej1...(yh)ejh, is applicable to the environment ejif it contains object x. When such rule is applied, object xpasses from ejto ej1,...,e jhpossibly modified into objects y1,...,y hrespectively. At any moment t(1 ≤t≤T) for each object xin environment ej, if there exist communication rules whose left-hand side is (x)ej, then one of these rules will be applied. If more than one such a rule can be applied to an object, the system selects one randomly, according to their probability which is given by pr(t). For each j(1 ≤j≤m) there is just one further restriction, concerning the consistency of charges: in order to apply several rules of RΠjsimultaneously to the same membrane, all the rules must have the same electrical charge on their right-hand side. Following the properties verified by the probabilistic functions, rules in Rand REcanbeclassifiedintoblocks of rules, as showed in definitions 2, 3 and 4. Definition 2. The left and right-hand sides of the rules are defined as follows: (a) Given a rule r∈Rof the form u[v]α i→u[v]α iwhere 1≤i≤q, α, α∈{0,+,−} and u, v, u,v∈Γ∗: –Theleft-handsideofris LHS(r)=(i, α, u, v). The charge of LHS(r)is charge(LHS(r)) = α. The length of LHS(r)is |u|+|v|, what indicates the cooperation degree of the rule. –The right-hand side of ris RHS(r)=(i, α,u ,v). The charge of RHS(r)is charge(RHS(r)) = α. The length of RHS(r)is |u|+|v|. (b) Given a rule r∈REof the form (x)ej pr −−−→ (y1)ej1···(yh)ejh,thelefthand side of ris LHS(r)=(ej,x), and the right-hand side of ris RHS(r)=(ej1,y 1)···(ejh,y h). Definition 3. Rules from Rcan be classified into consistent blocks associated with (i, α, α,u,v)as follows: Bi,α,α,u,v ={r∈R:LHS(r)=(i, α, u, v)∧charge(RHS(r)) = α} Definition 4. Rules from REcan be classified into (consistent) blocks associated with (ej,x)as follows: Bej,x ={r∈RE:LHS(r)=(ej,x)}. Recall that, according to the semantics of our model, the sum of probabilities of all the rules belonging to the same block is always equal to 1; in particular, rules with probability equal to 1 form individual blocks. Note that rules that have exactly the same left-hand side (LHS) belongs to the same block, but rules with overlapping (but different) left-hand sides are classified into different blocks. The latter leads to object competition, what is a critical aspect to manage with the simulation algorithms. Definition 5. Two blocks Bi1,α1,α 1,u1,v1and Bi2,α2,α 2,u2,v2are mutually consistent with each other, if and only if (i1=i2∧α1=α2)⇒(α 1=α 2). 3TheDCBA The goal of the DCBA (Direct distribution based on Consistent Blocks Algorithm) [17] is to perform a proportional distribution of objects among competing blocks (with overlapping LHS), determining in this way the number of times that each rule in m j=1 RΠj∪REis applied. I.e. the algorithm simulates the computational steps of a PDP systems. Algorithm 3.1 describes the main loop of the DCBA. It follows the same general scheme as its predecessors, DNDP and BBB [18] where the simulation of a computing step is structured in two stages: The first stage (selection), selects which rules are to be applied (and how many times) on each environment. The second stage (execution), implements the effects of applying the previously selected rules, yielding the next configuration of the PDP system. Note that, although every Πjhas the same set of rules R, the probability functions may be different for each environment. See [17] for a more detailed explanation and examples of how to apply this algorithm. As shown in Algorithm 3.1, the selection stage consists of three phases: Phase 1 distributes objects to the blocks in a certain proportional way, Phase 2 assures the maximality by checking the maximal number of applications of each block, and Phase 3 translates block applications to rule applications by calculating random numbers using the multinomial distribution. Algorithm 3.1. DCBA MAIN PROCEDURE Require: A Population Dynamics P system of degree (q, m), T≥1 (time units), and A≥1(Accuracy). The initial configuration is called C0. 1: INITIALIZATION (Algorithm 3.2) 2: for t←1toTdo 3: Calculate probability functions fr,j (t)andp(t). 4: C t←Ct−1 5: SELECTION of rules: •PHASE 1: distribution (Algorithm 3.3) •PHASE 2: maximality (Algorithm 3.4) •PHASE 3: probabilities (Algorithm 3.5) 6: EXECUTION of rules. (Algorithm 3.6) 7: Ct←C t 8: end for The INITIALIZATION procedure (Alg. 3.2) constructs a static distribution table Tjfor each environment. Two variables, Bj sel and Rj sel, are also initialized, in order to store the selected multisets of blocks and rules, respectively. Observation 1. Each column label of the tables Tjcontains the information of the corresponding block left-hand side. Observation 2. Each row of the tables Tjcontains the information related to the object competitions: for a given object, its row indicates which blocks are competing for it (those columns having non-null values). Algorithm 3.2. INITIALIZATION 1: Construction of the static distribution table T: –Column labels: consistent blocks Bi,α,α,u,v of rules from R. –Row labels: pairs (o, i)and(x, 0), for all object o∈Γ,x∈Σand membrane i, being 0 the identifier of the environments of the P system. –For each row labelled by (o, i) and column labelled by block Bi,α,α,u,v:place 1 kif oappears within i(0 ≤i≤q) with multiplicity kin the LHS of Bi,α,α,u,v. 2: for j=1to mdo (Construct the static expanded tables Tj) 3: Tj←T.(Initialize the table with the original T) 4: For each rule block Bej,x from RE, add a column labelled by Bej,x to the table Tj;placethevalue1atrow(x, 0) for that column. 5: Initialize the multisets Bj sel ←∅and Rj sel ←∅ 6: end for The distribution of objects among the blocks with overlapping LHS (competing blocks) is performed in selection Phase 1 (Algorithm 3.3). The expanded static tables Tjare used for this purpose in each environment, together with three different filter procedures. Filter 1 discards the columns of the table corresponding to non-applicable blocks due to mismatch charges in the LHS and in the configuration C t. Then, Filter 2 discards the columns with objects in the LHS not appearing in C t. Finally, in order to save space in the table, Filter 3discards empty rows. These three filters are applied at the beginning of Phase 1, and the result is a dynamic table Tt j(for the environment jand time step t). The semantics of the modeling framework requires a set of mutually consistent blocks before distributing objects to the blocks. For this reason, after applying Filters 1 and 2, the mutually consistency is checked. Note that this checking can be easily implemented by a loop over the blocks. If it fails, meaning that an inconsistency was encountered, the simulation process is halted, providing a warning message to the user. Nevertheless, it could be interesting to find a way to continue the execution by non-deterministically constructing a subset of mutually consistent blocks. Since this method can be exponentially expensive in time, it is optional for the user whether to activate it or not. Once the columns of the dynamic table Tt jrepresent a set of mutually consistent blocks, the distribution process starts. This is carried out by creating atemporalcopyofTt j, called TV t j, which stores the following products: –The normalized value with respect to the row: this is a way to proportionally distribute the corresponding object along the blocks. Since it depends on the multiplicities in the LHS of the blocks, the distribution, in fact, penalizes the blocks requiring more copies of the same object. This is inspired in the amount of energy required to gather individuals from the same species. –The value in the dynamic table (i.e. 1 k): this indicates the number of possible applications of the block with the corresponding object. –The multiplicity of the object in the configuration C t: this performs the distribution of the number of copies of the object along the blocks. Algorithm 3.3. SELECTION PHASE 1: DISTRIBUTION 1: for j=1to mdo (For each environment ej) 2: Apply filters to table Tjusing C t, obtaining Tt j. The filters are applied as follows: a. Tt j←T j b. Filter 1 (Tt j,C t). c. Filter 2 (Tt j,C t). d. Check mutual consistency for the blocks remaining in Tt j.If there is at least one inconsistency then report the information about the error, and optionally halt the execution (in case of not activating step 3.) e. Filter 3 (Tt j,C t). 3: (OPTIONAL) Generate a set St jof sub-tables from Tt j,formedbysetsof mutually consistent blocks, in a maximal way in Tt j(by the inclusion relationship). Replace Tt jwith a randomly selected table from St j. 4: a←1 5: repeat 6: for all rows Xin Tt jdo 7: RowSumX,t,j ←total sum of the non-null values in the row X. 8: end for 9: TV t j←Tt j(A temporal copy of the dynamic table) 10: for all non-null positions (X, Y )inTt jdo 11: multX,t,j ←multiplicity in C tat ejoftheobjectatrowX. 12: TV t j(X, Y )←multX,t,j ·(Tt j(X,Y ))2 RowSumX,t,j  13: end for 14: for all not filtered column, labelled by block B,inTt jdo 15: Na B←minX∈rows(Tt j)(TV t j(X, B)) (The minimum of the column) 16: Bj sel ←Bj sel +{BNa B}(Accumulate the value to the total) 17: C t←C t−LHS(B)·Na B(Delete the LHS of the block.) 18: end for 19: Filter 2 (Tt j,C t) 20: Filter 3 (Tt j,C t) 21: a←a+1 22: until (a>A)∨(all the selected minimums at step 15 are 0) 23: end for Algorithm 3.4. SELECTION PHASE 2: MAXIMALITY 1: for j=1to mdo (For each environment ej) 2: Set a random order to the blocks remaining in the last updated table Tt j. 3: for all block B, following the previous random order do 4: NB←number of possible applications of Bin C t. 5: Bj sel ←Bj sel +{BNB}(Accumulate the value to the total) 6: C t←C t−LHS(B)·NB(Delete the LHS of block B,NBtimes.) 7: end for 8: end for After the object distribution process, the number of applications for each block is calculated by selecting the minimum value in each column. This number is then used for consuming the LHS from the configuration. However, this application could be not maximal. The distribution process can eventually deliver objects to blocks that are restricted by other objects. As this situation may occur frequently, the distribution and the configuration update process is performed A times, where Ais an input parameter referring to accuracy. The more the process is repeated, the more accurate the distribution becomes, but the performance of the simulation decreases. We have experimentally checked that A=2givesthe best accuracy/performance ratio. In order to efficiently repeat the loop for A, and also before going to the next phase (maximality), it is interesting to apply Filters 2 and 3 again. After phase 1, it may be the case that some blocks are still applicable to the remaining objects. This may be caused by a low Avalue or by rounding artifacts in the distribution process. Due to the requirements of P systems semantics, a maximality phase is now applied (Algorithm 3.4). Following a random order, a maximal number of applications is calculated for each block still applicable. After the application of phases 1 and 2, a maximal multiset of selected (mutually consistent) blocks has been computed. The output of the selection stage has to be, however, a maximal multiset of selected rules. Hence, Phase 3 (Algorithm 3.5) passes from blocks to rules, by applying the corresponding probabilities (at the local level of blocks). The rules belonging to a block are selected according to a multinomial distribution M(N,g1,...,g l), where Nis the number of applications of the block, and g1,...,g lare the probabilities associated with the rules r1,...,r lwithin the block, respectively. Algorithm 3.5. SELECTION PHASE 3: PROBABILITY 1: for j=1to mdo (For each environment ej) 2: for all block BNB∈Bj sel do 3: Calculate {n1,...,n l}, a random multinomial M(NB,g 1,...,g l)with respect to the probabilities of the rules r1,...,r lwithin the block. 4: for k=1to ldo 5: Rj sel ←Rj sel +{rnk k}. 6: end for 7: end for 8: Delete the multiset of selected blocks Bj sel ←∅.(Useful for the next step) 9: end for Finally, the execution stage (Algorithm 3.6) is applied. This stage consists on adding the RHS of the previously selected multiset of rules, as the objects present on the LHS of these rules have already been consumed. Moreover, the indicated membrane charge is set. Algorithm 3.6. EXECUTION 1: for j=1to mdo (For each environment ej) 2: for all rule rn∈Rj sel do (Apply the RHS of selected rules) 3: C t←C t+n·RHS(r) 4: Update the electrical charges of C tfrom RHS(r). 5: end for 6: Delete the multiset of selected rules Rj sel ←∅.(Useful for the next step) 7: end for 4 The CUDA Programming Model With the commercial sector’s demands for video and gaming, it was foreseen by Elster [10] and others that graphics processor development would lead to devices suitable for High Performance Computing (HPC). With the introduction of NVIDIA’s CUDA [13] and AMD’s Stream SDK environments in 2007, the GPUs became more easily programmable and the era of GPGPU (General-Purpose Computing on Graphics Processing Units) truly began. GPUs can now be considered affordable computing solutions for speeding up computationally demanding applications. GPUs today typically have several hundred computational cores. By parallelizing and optimizing our codes for these cores, we have shown that GPUs can speed up applications ranging from smoothed hydrodynamics (SPH) [14] and real-time gradient vector flows to linear programming [21] and can even be used as an accelerator to compress I/O data for faster I/O speeds [1]. As was mentioned earlier, we have also successfully used GPUs to implement other P systems simulators [2,5,6]. 4.1 Compute Unified Device Architecture NVIDIA’s CUDA (Compute Unified Device Architecture) provides developers with a high-level programming model that allows developers to take full advantage of the NVIDIA’s powerful GPU hardware. To a CUDA programmer [23], the computing system is heterogeneous, consisting of a host (the CPU), and one or more massively parallel many-core devices (GPU systems). In modern software applications, there are often program sections that exhibit rich amount of data parallelism, a property where many arithmetic operations can be safely performed simultaneously on the program data structures. GPUs accelerate the execution of these applications by harvesting a large amount of data parallelism. CUDA is built around a scalable array of multithreaded Streaming Multiprocessors (SMs). The SM creates, manages, and executes concurrent threads (extra light weight processes) in hardware with with virtually no overhead [23]. This allows very fine-grained decomposition of problems by assigning, for instance, one thread to each data element. Each threads has a unique threadIdx, and is grouped in 1D, 2D, or 3D blocks which are organized in 1D or 2D grids. All the threads execute the same code, called kernel. Since all threads of a block are expected to reside on the same processor core and must share the limited memory resources of that core, current GPUs, The second test analyses the performance when increasing the parallelism level of the CUDA threads within thread blocks, that is, the number of rule blocks. The speedup achieved by pdp-gpu-sim versus pdp-omp-sim is showed in Figure 4. The number of simulations is fixed to 50, and the environments to 20 (hence, a total of 1000 thread blocks). The number of objects is proportionally increased together with the number of rule blocks, in such a way that the ratio for number of rule blocks and number of objects is always 2. The mean LHS length is 1.5 (this is typical for many real ecosystem models, as seen in the literature). The speedup gets stable to around 7x on the number of rule blocks for the GPU versus CPU. For the multicore versions with 2 and 4 CPUs, the speedups are maintained to 4.3x and 3x, respectively. In our experiments, this number is also achieved when running with 1000000 rule blocks.                                    Fig. 4. Scalability of the simulators when increasing the number of rule blocks The third test is for the second parallelism level in CUDA, concerning thread blocks. It is directly related with the number of environments and simulations. The result is shown in Figure 5. In this experiment, the number of rule blocks is fixed to 10000, the number of objects to 7024 and the mean LHS length is 2. The number of environments is fixed to 1 when increasing the simulations, and vice versa. As it can be seen, for low values, the speedup is demoted below 1. These values come from the fact of insufficient number of thread blocks to fulfill the GPU resources. Another trend shown is that as the number of simulations increases, the advantage of parallelizing by simulations increases. The same effect is observed for environments. This trend is stabilized to 3.5x for high values. However the parallelism over simulations is better carried out by the GPU, giving lower speedups for environments. As stated in [16], parallelizing by simulations yields the largest speedups on multicore platforms. Therefore, we finalize the first benchmark by comparing these results with the GPU. Rule blocks are fixed to 50000, environments to                                          Fig. 5. Scalability – increasing the number of simulations and environments                                   ! Fig. 6. Scalability of the simulators when increasing the number of simulations 20, objects to 5000 and mean LHS length to 1.5. As shown in Figure 6, the GPU achieves better runtime than the multicore implementations. The speedup is maintained to 4.5x using one core, 3.5x for 2 cores, and 2.7x for 4 cores. The results of the second benchmark are shown in Table 1. This profile has been calculated running the simulator with 10000 rule blocks, 20 environments, 50 simulations, 5000 objects and two different mean LHS lengths, 1.5 (test A) and 3 (test B). Phase 1 is the most complex part in the simulation (taking more than 50% of the runtime on the CPU). In test A, the GPU implementation offers for phase 1 up to 14x of speedup. Therefore, the percentage of the execution time is decreased to 30%. Following with Test A, Phase 2 takes only the 12% of the execution time on CPU. However, the GPU can only accelerate this phase by 2x. Therefore, this phase becomes the most expensive when executing the simulator on the GPU (47%). Our novel implementation, ph2-simorder-dyncomp, is close (timewise) to the sequential implementation. Indeed, as mentioned above, this phase Table 1. Profiling the simulators for GPU and 1 core CPU Test A (mean LHS length 1.5) Test B (mean LHS length 3) %CPU %GPU Speedup %CPU %GPU Speedup Phase 1 53.7% 30.1% 14.23x 55.3% 12% 8.52x Phase 2 12.6% 47% 2.13x 18.4% 82.8% 0.4x Phase 3 22.6% 13.7% 13.2x 14% 2.2% 11.72x Phase 4 11.1% 9.2% 9.7x 12.3% 3% 7.43x is the most challenging to parallelize. Special efforts have to be considered here. On the other side, Phases 3 and 4 are relatively lightweight, and are successfully accelerated (up to 13x and 9.7x, respectively). Hence, our library random binomial generation based on CuRAND is well suited for Phase 3. Finally, as shown in Figure 3, the performance of pdp-gpu-sim decreases as the mean LHS length is increased. For Test B, the overall speedup decreases from 7.9x (Test A) to 1.8x (Test B). The percentage of time consumed by Phase 2 is dramatically increased for the GPU, taking up to 83%. Thus, the competition degree of rule blocks is a limiting factor in performance, which fully correlates with the achieved results. 7 Conclusions and Future Work In this paper we have presented the first GPU-based version of a simulator for PDP systems with CUDA. We improved the memory utilization of both our GPU-based version and our previous OpenMP-based version that we benchmarked against. We benchmarked a set of randomly generated PDP systems (without biological meaning), achieving speedups of up to 7x for large sizes on NVIDIA Tesla C1060 GPU over our multi-core version. We used a general CUDA design for the GPU part of our simulator: environments and simulations are distributed through thread blocks, and rule blocks among threads. Phases 1, 3 and 4 were efficiently executed on the GPU, however Phase 2 was poorly accelerated, since it is inherently sequential. For future work, the simulator is going to be reconnected with the P-Lingua framework. The pLinguaCore library is able to parse P-Lingua files, simulate P systems computations, and translate P-Lingua files to other file formats. In this respect, a new file format is going to be designed to become the input of our simulator. That way, the same P-Lingua files can contain the input data for the Java simulator inside pLinguaCore and the CUDA simulator. Moreover, we would like to validate our simulator with real ecosystems models. This would also help us to adopt model-oriented heuristics to improve the CUDA design. Finally, it would be useful to design a communication protocol and a Graphics User Interface in order to connect the simulation pipeline with the end-user. This simulator represents our first attempt on simulating PDP systems using the GPU. Hence, improvements can be done to our new CUDA designs (i.e. iterating from objects to blocks (gather strategy)), optimizing the code, and/or making more efforts on Phase 2. Future work will also take advantage of the new GPU architectures, such as the NVIDIA Kepler, and AMD APU, which deliver better performance than the Tesla C1060 we used for our experiments here. Acknowledgments. The authors acknowledge the support of “Proyecto de Excelencia con Investigador de Reconocida Val´ıa” of the “Junta de Andaluc´ıa” under grant P08-TIC04200, and the support of the project TIN2009-13192 of the “Ministerio de Educaci´on y Ciencia” of Spain, both co-financed by FEDER funds. NVIDIA’s support and donations to our HPC-Lab at NTNU is also gratefully acknowledged. References 1. Aqrawi, A.A., Elster, A.C.: Bandwidth Reduction through Multithreaded Compression of Seismic Images. In: IEEE International Symposium on Parallel & Distributed Processing (IPDPSW 2011), pp. 1730–1739 (2011) 2. Cabarle, F., Adorna, H., Mart´ınez-del-Amor, M.A., P´erez-Jim´enez, M.J.: Improving GPU Simulations of Spiking Neural P Systems. Romanian Journal of Information Science and Technology 15(1), 5–20 (2012) 3. Cardona, M., Colomer, M.A., Margalida, A., Palau, A., P´erez-Hurtado, I., P´erez- Jim´enez, M.J., Sanuy, D.: A computational modeling for real ecosystems based on P systems. Natural Computing 10(1), 39–53 (2011) 4. Cardona, M., Colomer, M.A., Margalida, A., P´erez-Hurtado, I., P´erez-Jim´enez, M.J., Sanuy, D.: A P System Based Model of an Ecosystem of Some Scavenger Birds. In: P˘aun, G., P´erez-Jim´enez, M.J., Riscos-N´u˜nez, A., Rozenberg, G., Salomaa, A. (eds.) WMC 2009. LNCS, vol. 5957, pp. 182–195. Springer, Heidelberg (2010) 5. Cecilia, J.M., Garc´ıa, J.M., Guerrero, G.D., Mart´ınez-del-Amor, M.A., P´erez- Hurtado, I., P´erez-Jim´enez, M.J.: Simulation of P systems with Active Membranes on CUDA. Briefings in Bioinformatics 11(3), 313–322 (2010) 6. Cecilia, J.M., Garc´ıa, J.M., Guerrero, G.D., Mart´ınez-del-Amor, M.A., P´erez- Hurtado, I., P´erez-Jim´enez, M.J.: Simulating a P system based efficient solution to SAT by using GPUs. Journal of Logic and Algebraic Programming 79(6), 317–325 (2010) 7. Cheruku, S., P˘aun, A., Romero-Campero, F.J., P´erez-Jim´enez, M.J., Ibarra, O.H.: Simulating FAS-induced apoptosis by using P systems. Progress in Natural Science 17(4), 424–431 (2007) 8. Colomer, M.A., Lav´ın, S., Marco, I., Margalida, A., P´erez-Hurtado, I., P´erez- Jim´enez, M.J., Sanuy, D., Serrano, E., Valencia-Cabrera, L.: Modeling Population Growth of Pyrenean Chamois (Rupicapra p. pyrenaica) by Using P-Systems. In: Gheorghe, M., Hinze, T., P˘aun, G., Rozenberg, G., Salomaa, A. (eds.) CMC 2010. LNCS, vol. 6501, pp. 144–159. Springer, Heidelberg (2010) 9. Colomer, M.A., P´erez-Hurtado, I., P´erez-Jim´enez, M.J., Riscos, A.: Comparing simulation algorithms for multienvironment probabilistic P system over a standard virtual ecosystem. Natural Computing, doi:10.1007/s11047-011-9289-2 10. Elster, A.C.: High-Performance Computing: Past, Present, and Future. In: Fagerholm, J., Haataja, J., J¨arvinen, J., Lyly, M., R˚aback, P., Savolainen, V. (eds.) PARA 2002. LNCS, vol. 2367, pp. 433–444. Springer, Heidelberg (2002) 11. Garc´ıa-Quismondo, M., Guti´errez-Escudero, R., P´erez-Hurtado, I., P´erez-Jim´enez, M.J., Riscos-N´u˜nez, A.: An Overview of P-Lingua 2.0. In: P˘aun, G., P´erez-Jim´enez, M.J., Riscos-N´u˜nez, A., Rozenberg, G., Salomaa, A. (eds.) WMC 2009. LNCS, vol. 5957, pp. 264–288. Springer, Heidelberg (2010) 12. Kachitvichyanukul, V., Schmeiser, B.W.: Binomial random variate generation. Communications of the ACM 31(2), 216–222 (1988) 13. Kirk, D., Hwu, W.: Programming Massively Parallel Processors: A Hands on Approach, MA, USA (2010) 14. Krog, Ø.E., Elster, A.C.: Fast GPU-Based Fluid Simulations Using SPH. In: J´onasson, K. (ed.) PARA 2010, Part II. LNCS, vol. 7134, pp. 98–109. Springer, Heidelberg (2012) 15. Harris, M.: Mapping computational concepts to GPUs. In: ACM SIGGRAPH 2005 Courses, NY, USA (2005) 16. Mart´ınez-del-Amor, M.A., Karlin, I., Jensen, R.E., P´erez-Jim´enez, M.J., Elster, A.C.: Parallel Simulation of Probabilistic P Systems on Multicore Platforms. In: Proceedings of the Tenth Brainstorming Week on Membrane Computing, vol. II, pp. 17–26 (2012) 17. Mart´ınez-del-Amor, M.A., P´erez-Hurtado, I., Garc´ıa-Quismondo, M., Mac´ıas- Ramos, L.F., Valencia-Cabrera, L., Romero-Jim´enez, A., Graciani, C., Riscos- N´u˜nez, A., Colomer, M.A., P´erez-Jim´enez, M.J.: DCBA: Simulating Population Dynamics P Systems with Proportional Object Distribution. In: Proceedings of the Tenth Brainstorming Week on Membrane Computing, vol. II, pp. 27–56 (2012) 18. Mart´ınez-del-Amor, M.A., P´erez-Hurtado, I., P´erez-Jim´enez, M.J., Riscos-N´u˜nez, A., Sancho-Caparrini, F.: A simulation algorithm for multienvironment probabilistic P systems: A formal verification. International Journal of Foundations of Computer Science 22(1), 107–118 (2011) 19. P˘aun, G.: Computing with membranes. Journal of Computer and System Sciences 61(1), 108–143 (2000), TUCS Report No 208 20. P˘aun, G., Rozenberg, G., Salomaa, A. (eds.): The Oxford Handbook of Membrane Computing. Oxford University Press (2010) 21. Spampinato, D.G., Elster, A.C.: Linear optimization on modern GPUs. In: IEEE International Symposium on Parallel & Distributed Processing (IPDPS 2009), pp. 1–8 (2009) 22. Terrazas, G., Krasnogor, N., Gheorghe, M., Bernardini, F., Diggle, S., C´amara, M.: An Environment Aware P-System Model of Quorum Sensing. In: Cooper, S.B., L¨owe, B., Torenvliet, L. (eds.) CiE 2005. LNCS, vol. 3526, pp. 479–485. Springer, Heidelberg (2005) 23. NVIDIA CUDA programming guide 4.0 (2012), http://www.nvidia.com/cuda 24. The P-Lingua web page, http://www.p-lingua.org