Efficient simulation execution of cellular automata on GPU

Simulation Modelling Practice and Theory 118 (2022) 102519 Available online 5 March 2022 1569-190X/© 2022 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Contents lists available at ScienceDirect Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat Efficient simulation execution of cellular automata on GPU Daniel Cagigas-Muñiz∗, Fernando Diaz-del-Rio, Jose Luis Sevillano-Ramos, Jose-Luis Guisado-Lizar Department of Computer Architecture and Technology, Universidad de Sevilla, Avenida Reina Mercedes s/n, 41012 Sevilla, Spain ARTICLE INFO Keywords: Cellular automata Parallel computing Performance optimization Stencil computation Graphics Processing Units ABSTRACT Graphics Processing Units (GPUs) can be used as convenient hardware accelerators to speed up Cellular Automata (CA) simulations, which are employed in many scientific areas. However, an important set of CA have performance constraints due to GPU memory bandwidth. Few studies have fully explored how CA implementations can take advantage of modern GPU architectures, mainly in the case of intensive memory usage. In this paper, we make a thorough study of techniques (stencil computing framework, look-up tables, and packet coding) to efficiently implement CA on GPU, taking into account its detailed architecture. Exhaustive experiments to validate these implementation techniques for a number of significant memorybounded CA are performed. The CA analysed include the classical Game of Life, a Forest Fire model, a Cyclic cellular automaton, and the WireWorld CA. The experimental results show that implementations using the presented techniques can significantly outperform a baseline standard GPU implementation. The best performance results of all known implementations of memory bounded CA were obtained. Moreover, some of the techniques, like look-up tables or temporal blocking, are indeed relatively easy to implement or to apply when the transition rules are simple. Finally, detailed descriptions and discussions of the indicated techniques are included, which may be useful to practitioners interested in developing high performance simulations in efficient languages based on CA on GPU. 1. Introduction Cellular Automata (CA) are mathematical models that evolve in discrete time steps, composed of elements called ‘‘cells’’ that change their state at each time step according to some rules, by taking into account the state of the neighbouring cells. CA are commonly used in complex and/or dynamic system modelling and simulation, mainly in the fields of physics, biology, and computer science [1–4]. A sequential implementation of a cellular automaton using slow interpreted languages can be enough in certain situations in which the output time response is not critical, the CA model evolution rules are simple, the CA data are relatively small, and/or experiments do not need to be repeated many times. In these cases, many practitioners just prototype for slow CPU performance nonparallel programming languages like Python, Matlab, or Octave. However, there are many real-world applications of CA simulation models that demand a high level of computational capability because their CA model is complicated in terms of evolution rules or data usage, for real-time (or time-constrained) applications, or for parametric studies in which the same experiment must be executed many times under different initial configurations (in ∗Corresponding author. E-mail addresses: [email protected] (D. Cagigas-Muñiz), [email protected] (F. Diaz-del-Rio), [email protected] (J.L. Sevillano-Ramos), [email protected] (J.-L. Guisado-Lizar). https://doi.org/10.1016/j.simpat.2022.102519 Received 8 January 2022; Received in revised form 22 February 2022; Accepted 24 February 2022 Simulation Modelling Practice and Theory 118 (2022) 102519 2 D. Cagigas-Muñiz et al. which case it is critical to accelerate every individual experiment). Therefore, these applications of CA require the usage of parallel computing techniques. Examples involve scientific problems as diverse as fluid dynamics [5], cancer growth [6], 3D metal forming [7], solidification in materials science [8], or urban land use simulations [9], among many others. It is therefore interesting to analyse some programming fundamentals on how CA implementations on modern computing platforms can be improved when performance is important. Computer-based implementations of CA are inherently parallel because each cell can be processed independently at each time step. With the arrival in 2007 of programming languages for Graphical Processing Units (GPUs) such as Cuda and later OpenCL, the design and implementation of CA experienced a new impetus. The internal architecture of GPUs has many more computational units (processors) than a CPU. This feature makes GPUs ideal hardware accelerators for implementing CA. There is much room for performance improvement, especially for memory-bounded CA, such as the GoL. Memory-bounded CA are cellular automata that use more time in memory accesses than in computing cell states for the next time step. This means that the GPU arithmetic computational units are not being efficiently used. At each time step of CA, cell computations can be performed in parallel by several processors without the need for communication or synchronization among them. Several experimental works can be found in the scientific literature about the performance benefits of using GPUs over CPUs (see Section 2). There is, however, one problem when coding performance-sensitive algorithms using GPUs. GPU manufactures do not provide access to the machine code. For example, in the case of NVIDIA Cuda programming language, only Parallel Thread Execution (PTX) is provided. PTX is a low-level parallel thread execution virtual machine and instruction set architecture (ISA) [10]. Although the syntax of PTX is similar to that of an assembly level programming language, there is not a one to one correspondence between a PTX instruction and a GPU machine code instruction. Moreover, GPU’s instruction set architecture (ISA) and/or instruction encoding changes from architecture to architecture. Another important related problem are compilers. In the case of NVIDIA and Cuda, the compiler is a ‘black box’ that applies code optimization. Nevertheless, it is not possible to determine whether compiler optimizations are always appropriate. Therefore, there is no way of knowing exactly how the high-level code (Cuda or OpenCL) is actually executed on the GPU. Only some general guidelines based on the GPU architecture are indicated by manufacturers. This makes experimentation, profiling, and testing necessary to ensure effective code improvement. CA are not an exception despite its apparent simplicity. Getting the maximum performance from GPU implementations is not a trivial task. Knowledge on low level hardware architecture details is needed to take advantage of modern GPU full capabilities. Even GPU programmers do not always squeeze the full potential of their hardware. Concurrent and/or parallel programming has obvious advantages, but it is always more complex to implement. However, as mentioned before, high performance is desired. For a better comprehension of the main bottleneck involved in GPU performance, our analysis is clearly extensible to higher dimensional problems. Nonetheless, many code optimization techniques described in this study can be extended to three-dimensional CA more or less directly. Three-dimensional CA are another good example of massive data structures that need an efficient GPU implementation for a better performance. The contributions of this paper can be summarized in three points: 1. An analysis of which GPU architecture aspects influence most in the performance of CA implementations. 2. An analysis of different new algorithms and/or techniques to improve baseline CA implementations performance on GPUs, breaking the current GPU performance barrier for some class of CA. 3. The implementation and validation (with experimental results) of the new algorithms and/or techniques proposed. As a result, a new highly efficient implementation technique is proposed for memory bandwidth bounded CA. This article is organized as follows. Related works are discussed in Section 2. A generic CA baseline implementation on GPU is described in Section 3. Performance issues related to the implementation of CA on GPU architectures are analysed in Section 4. Section 5presents new coding techniques for higher performance, which are based on previous analysis. The proposed techniques are tested using some CA models in Section 6. The experimental results are discussed in Section 7. Finally, conclusions and future work are outlined in Section 8. 2. Related works There have been some studies on the application of GPUs for CA in the scientific literature. Most of them are focused on highlighting the benefits of computing CA on GPUs versus CPUs. Examples can be found in [11,12], and [13]. In these studies, there is little discussion about GPU optimization algorithms and/or techniques applicable to CA. In addition, they usually focus on only one cellular automaton: the Game of Life (GoL). GoL was created by Conway [14] in 1970 and it is the most referenced and well-known cellular automaton. Unfortunately, these references did not provide data on more (complex) CA models. GoL is considered a ’toy model’ cellular automaton that is often far from real dynamic models. An exception to this practice can be found, for example, in [15] where laser dynamics is modelled by using a cellular automaton employing high performance multiprocessors and GPUs. Several algorithms and issues that influence the CA performance in GPUs were studied in [16] using again GoL cellular automaton as an example. Issues such as the influence of shared memory on GPUs and the thread block size on Cuda were analysed. In particular, this work did not find any influence of the block size on the algorithms used; thus these authors set always the block size to 32 ×32 threads. They also tested an algorithm called ‘‘multicell algorithm’’, in which each GPU thread processes two cells. Although it can Simulation Modelling Practice and Theory 118 (2022) 102519 3 D. Cagigas-Muñiz et al. Fig. 1. Pseudo code description of a generic cellular automaton program in GPU. be beneficial in other types of algorithms or models, it did not improve the results of the baseline GoL implementation. This study provided several interesting results but left significant room for further research to improve the CA performance on modern GPUs. No solutions to improve the baseline GoL implementation using GPUs were provided. Previously, there were also some nonspecific studies in [17]. However, some of the results and conclusions provided were totally opposed to [16]. This is due to the fact that the authors in [17] worked with old NVIDIA graphic cards (Fermi architecture). These graphics cards were the first generation of full Cuda programmable GPUs that did not have any kind of cache memory in their internal architecture. As it will be seen in later sections, this issue is absolutely critical to understand the performance of CA on modern GPUs. 3. Cellular automata model The typical bi-dimensional CA (Cellular Automata) workflow on a GPU consists of first declaring two grids (arrays) 𝐴and 𝐵in memory. The first grid 𝐴represents the current state of the CA and the second grid 𝐵the state after a time or computation step. Note that using only one grid for the current and next state would produce incorrect values for GPU execution, since a thread could modify a cell state before the other thread has read this state. The content of 𝐴(current state) is prefilled with the initial data of the CA on the host computer (usually a multiprocessor). Then, its content is transferred to the GPU memory. From that moment on, both grids (arrays) are only accessed in the GPU memory. In each time iteration or computation step, the grid 𝐵(next state) is calculated based on the content of matrix 𝐴(current state). At the end of the computation step, the roles of 𝐴and 𝐵are exchanged to prepare the next time iteration. In programming languages like C/C++, this implies the use of pointers (𝑝_𝑐𝑢𝑟𝑟𝑒𝑛𝑡_𝑠𝑡𝑎𝑡𝑒,𝑝_𝑛𝑒𝑥𝑡_𝑠𝑡𝑎𝑡𝑒). Cuda and OpenCL are based on C/C++ because they are programming languages that focus on code efficiency. Once the last iteration or computation step has been completed, the final contents of the grid containing the current state of the CA must be transferred from the GPU (device) memory to the CPU (host) memory. Fig. 1 shows a pseudo-algorithm that summarizes a typical CA implementation on a GPU. The 𝑐𝑜𝑚𝑝𝑢𝑡𝑒_𝑠𝑡𝑒𝑝 function from Fig. 1 implies the execution of at least one kernel on the GPU. A kernel is nothing more than a function executed within the GPU. One of the critical aspects of obtaining good performance is minimizing transfers between the GPU and the host computer as they involve memory latency. This is why memory transfers between GPU and CPU are usually limited to the beginning and the end of the CA simulation. For the same reason, the number of iterations or computation steps performed should be relatively high to take advantage of the parallel architecture of the GPU. The use of relatively small 𝐴and 𝐵 grids of data, the continuous transfer of information between the GPU and CPU, and/or a few time steps can result in an important decrease of performance. 4. Analysis of CA performance As stated in the Introduction, there is much room for time enhancement in memory-bounded CA. A preliminary analysis was done by us in [18], where we show that the GoL implementations are far from taking full advantage of current architectures using the roofline model [19]. This model is an excellent tool for detecting promptly potential bottlenecks and performance issues. Going forward, in this section, the two factors that have a major influence on the CA performance are to be analysed in more depth. Simulation Modelling Practice and Theory 118 (2022) 102519 4 D. Cagigas-Muñiz et al. Fig. 2. Comparison of performance of GPU memory types (global, shared and texture memories) in GoL cellular automaton and different grid sizes. A NVIDIA GTX 1650 TI GPU (Turing architecture) was used. 4.1. GPU memory model management The objective in this subsection is to determine which GPU memory management is more effective for CA implementation. GPUs have several types of memory and it is responsibility of the programmer to determine its management. This first analysis will provide a first approach of how CA should be coded. To do this, several CUDA implementations of GoL (Game of Life) cellular automaton were tested. There are only two states in GoL: dead or alive. Alive cells remain alive if there are two or three alive neighbour cells. Dead cells become alive with three alive neighbour cells. In any other case, cells become dead or remain dead. Moore neighbourhood (up, down, left, right, and four diagonals) is taken into account. There are some variations of this cellular automaton, including even 3D versions. The GoL is the base for much more complex CA. Cuda is the GPU programming language used in every experiment and implementation of this paper. It is the native programming language of the NVIDIA manufacturer and it is the predominant one currently. The concepts explained here are perfectly valid and analogous to other GPU languages such as OpenCL. From the programmer’s point of view, cuda allows the use of three main types of memory: registers (local memory), shared memory, and global memory. The general strategy for achieving a good acceleration/performance in GPUs is to try to avoid access to global memory as much as possible. Copying data from global memory to shared memory or using textures are alternatives to avoid many global memory accesses. Textures are an additional type of memory present in NVIDIA GPUs that are specifically geared towards image processing. Thus, three implementations of the GoL based on [20] were tested: using global memory, using shared memory, and using textures. To study the effect of memory, these three memory models have also been tested for different CA (grid) sizes. The summary of the results obtained are shown in Fig. 2 for a total of 1000 time steps. The overall conclusion that can be drawn from these first results is that although textures are ideal for processing 2-dimensional images, they are not so ideal for implementing bi-dimensional CA. As for the implementation of the GoL using global vs. shared memory at the thread block level, theory indicates that using shared memory can reduce access to global memory if data sharing among threads is carefully programmed. The global memory data bus acts as a bottleneck in case all threads want to access it. However, it is interesting to remark that GoL experimental results in Fig. 2 do not confirm this statement. Results using only global memory are slightly worse than those of the shared memory case. This effect is more noticeable for larger (grid) CA. A detailed analysis of the CA model (see previous Section 3) and the GPU architecture itself can provide the key to explain this behaviour. The cellular automaton grid (data matrix) that has the current state of the CA only needs to be read from the global memory for each time/computation step. In addition, this data structure is not modified during the entire computation step. This data, like the automatic variables defined in a thread block, are stored in cache memory and in registers within a Multiprocessor Streaming (SM). Moving data that is read from global memory to shared memory in a block of threads creates an unnecessary overhead because that information is already cached either in registers or in the 1 or 2 cache memory level. Write memory access has a similar behaviour. The cell information processed by a block of threads is written only once and in a different area of the global memory. Therefore, moving the information from global memory to shared memory in each block of threads, and then again to global memory, does not reduce the data traffic to global memory. It only generates a greater overhead by having to move data unnecessarily from one type of memory to another. Cache memory is again a key component because it groups writes in blocks without programmer intervention. Shared memory is useful to gain performance only when threads in a block have to share data frequently over the same memory area (i.e., read/write several times to the same data in the same time step), before writing the results to global memory. This is not the case of the GoL and many CA. Simulation Modelling Practice and Theory 118 (2022) 102519 5 D. Cagigas-Muñiz et al. Fig. 3. Acceleration obtained when changing the GoL cellular automaton cell data type from ‘‘int’’ to ‘‘char’’. Fifteen independent executions (tests) for each grid size and data type were performed. The mean time in milliseconds was taken. An NVIDIA GTX 1650 TI graphics card was used. Only the main GoL Cuda kernel was measured using the NVIDIA Profiler. The main Cuda kernel run time is between 74% (256 grid size) and 99% (16384 grid size) of the total program execution time on the GPU. As a first conclusion, it can be stated that, having enough internal registers and a good cache memory system (which is being improved with each new GPU generation), CA implementations tend to read and write each cell only once from the global memory in each time/computation step. This has also a positive aspect: programming using only global memory is easier, and therefore Cuda code is cleaner. The negative effect of using shared memory when programming CA in the absence of cache memory was also briefly commented in [16]. 4.2. Cell size Another factor to study involves reducing cell size. In the experiments of Section 4.1, the data type ‘‘int’’ was used for each cell. That means 32 bits per cell in memory. This allows to code up to 232 possible states per cellular automaton (grid) cell. In the case of the GoL, only two states per cell are possible (live or dead). Therefore, an 8-bit data type ‘‘char’’ per cell is enough for the GoL cellular automaton and many CA. This reduces memory usage by 4 times with a consequent decrease in memory latency. There are studies of GoL on multiprocessors such as [21] in which cells are encoded with a single bit. Nonetheless, apart from the technical complexity, this solution is not generalizable to other CA. Because Cuda (and OpenCL) is based on C/C++ programming language, it is not possible to use data types smaller than 8 bytes (‘‘char’’ type). The Fig. 3 shows the acceleration factor obtained by changing the cell type from ‘‘int’’ to ‘‘char’’ in the GoL cellular automaton. It can be clearly seen that: 1. The acceleration is very irregular for small grid sizes: 1% for a 256 ×256 grid size and 41% for a 512 ×512 grid size. 2. Then, the performance improvement tends to be constant for large grids. This is determined by memory bandwidth constraints and GPU functional units. The maximum acceleration obtained is low in proportion to the decrease in memory usage. GPUs are optimized to work with 32-bit floating point data, which is the same size as the ‘‘int’’ data. Although the 8-bit ‘‘char’’ data type represents a decrease of x4 data in memory, the maximum performance obtained is below 20% for large grid sizes (see Fig. 3). GPUs have functional units designed for float or double data types that involve 32 or 64 bits respectively. The 8 bit char data type is not well suited for the GPU functional units, and can have negative effects in relatively small CA grids as it can be observed in Fig. 3. Nevertheless, in the case of large CA, codifying each CA cell with a smaller data size variable reduces the overall data size and therefore the performance should be better. This gives a hint to readdress memory accesses by using standard and wide variables, as it will be presented in Section 5.2. 5. Techniques for runtime-efficient coding of CA Based on the previous works mentioned in Section 2, CA model description of Section 3, and the analysis results of Section 4, some techniques for promoting CA efficient performance are proposed. Some of these techniques have been used before, but no study has been carried out to extend them to generic CA on GPUs. Simulation Modelling Practice and Theory 118 (2022) 102519 6 D. Cagigas-Muñiz et al. 5.1. Look-up tables An interesting performance improvement consists of trying to reduce or eliminate selective structures (𝑖𝑓 −𝑡ℎ𝑒𝑛 −𝑒𝑙𝑠𝑒). GPUs execute the same instructions for all threads (basic computational units) in a warp or block. This style of programming, called single-instruction multiple-thread (SIMT), works well when all threads take the 𝑖𝑓 −𝑡ℎ𝑒𝑛 path at once. Otherwise, several passes of the block/warp execution are required. One pass will be needed for those threads that follow the 𝑡ℎ𝑒𝑛 part and another pass for those that follow the 𝑒𝑙𝑠𝑒 part. These passes are sequential to each other, thus incrementing the execution time [22]. In general, the 𝑖𝑓 −𝑡ℎ𝑒𝑛 −𝑒𝑙𝑠𝑒 or 𝑠𝑤𝑖𝑡𝑐ℎ −𝑐𝑎𝑠𝑒 selective structures appear naturally when coding CA rules and it is convenient to suppress them. Programmers usually do not take this issue into account when coding CA. An alternative is to use what it is going to be defined in this study as ‘‘Look-Up Tables’’ (LUTs) [18]. These tables are data structures that have as input the current state of a cell (row coordinates) and the state and/or number of neighbouring cells (column coordinates). They are coded as bi-dimensional arrays. For instance, in the case of the GoL, this involves using a matrix of 2 rows by 9 columns. The rows represent the two possible states of a cell (live or dead) and the columns represent the number of neighbouring live cells. The content of this matrix indicates the state at which the cell should change in the next time/computation step. LUTs coding and management is very dependent on the cellular automaton. Nevertheless, in the case of CA with intensive memory usage, the transition rules between states tend to be relatively simple (i.e., need low computational resources). LUTs can be simplified in this set of CA and are relative easy to code. This issue will be analysed in detail with practical CA in Section 6. 5.2. Packet coding A plausible solution to the limitation imposed by the GPU bandwidth is to access memory only in 𝑓𝑙𝑜𝑎𝑡 or 𝑑𝑜𝑢𝑏𝑙𝑒 data chunks. Therefore, each thread computes more than one cell. In this way, the computational load of each thread is increased. As mentioned in Section 2, a somewhat related but different approach was analysed and named ‘‘multicell algorithm’’ in [16], in which each thread (cell) accesses to its neighbours sequentially (see Fig. 4, middle). However, memory accesses were done without taking into account the GPU memory data bus architecture. As a consequence, these authors obtained slightly worse results even than the GoL Baseline version. Although the original idea was correct, because the ratio of read accesses per cell is reduced from 9/1 (3 ×3 for a cell) to 12/2 (4 ×3 for two contiguous cells), the GPU architecture is still not really used efficiently. When a thread processes two cells of a CA, it produces two different (non-coalesced) access requests to the memory (see Fig. 4, middle). In contrast, it is possible to read/write a set of cells packed in just one memory access (see Fig. 4, bottom). Therefore, memory access requests can be reduced by x2, x4, or more. To develop this algorithm, several cells must be codified in a ‘‘supercell’’ whose size is that of the Functional Unit bus width. In the case of programming languages like C/C++ (and cuda), the maximum data size available is 64 bits. The ideal scenario, however, would be to utilize even higher data sizes. Fig. 5 shows how to code 8 cells of the same type in a supercell for the case of GoL computation. It can be observed that, although 8 neighbouring supercells must be read to calculate the new states of a particular supercell, central cells (that is, all except the first and the last one) only need to access cells from the upper and lower supercells. This technique, which is going to be called ‘‘Packet Coding’’ in this paper, has been scarcely used and for very specific cases. To the authors’ knowledge, the unique deep study in the field of CA is [21] and only multiprocessors were used. In [21] the GoL was coded with 1 bit per cell for an CPU implementation. The performance results were remarkable; however, as mentioned in Section 4.2 this solution is not extensible to every CA. Experiments using Packet Coding for different CAs are described in Section 6. 5.3. Temporal blocking A third approach for improving CA performance on GPUs is to take advantage of the studies that have been made in stencil computing [23]. Stencil computation has been used in many applications and ranges from simulation to machine vision, machine learning applications, or partial differential equation (PDE) solving. This has led to a great deal of research. A stencil code consists of updating the elements of an array (2D or 3D) based on a fixed pattern. These array elements are also often called cells. The similarities with CA are obvious and GPUs have also been used to achieve higher performance [24,25]. Part of the solution to the low CA performance is to improve the utilization of the GPU computer units. In stencil computing, the CA grid is usually divided into subplanes of the same size (named tiles in stencil computation) and each subplane or tile is processed by a thread block. At each time/computation step, all tiles of the grid are fully processed. Going further, another more complex technique for stencil computation is temporal blocking, which enhances the temporal reuse of data, and thus reduces the number of data transfers from/to the global GPU memory. It is based on a time-tiled execution, but the operations from several consecutive time steps are combined to exploit data reuse in memory [26]. However, the cell values that are in the edges (ghost or halo zones) depend on other tiles. In the case of GPUs, this means a dependency between thread blocks. This implies that in temporal blocking there must be periodic synchronization between threads. In Fig. 6 a basic scheme of how temporal blocking could be applied to a CA divided into 4 tiles is shown. This technique applied to 3-dimensional stencils is called 3.5D blocking in [27]. Global memory bandwidth reduction is thus achieved by combining several time/computation steps of the same tile to be executed consecutively. Intermediate data reside in registers, shared memory, or cached memory. The problem when using GPUs Simulation Modelling Practice and Theory 118 (2022) 102519 7 D. Cagigas-Muñiz et al. Fig. 4. Memory access patterns comparison for different CA codifications when accessing the North neighbours of a warp of cells (likewise for the rest of neighbours). From up to bottom: in a classical implementation looking for the north cell provokes a coalesced access to a set of consecutive bytes; on the contrary, the multicell coding introduces non-coalesced accessing (only one out of two elements in the case of a two-cell per thread implementation); finally, the novel packet coding reproduces again coalescing accesses but for wider element sizes. Fig. 5. Packet coding example: one supercell codes 8 cells, and 8 supercells are needed to compute a new supercell in a cellular automaton using Moore neighbourhood. is that, unlike CPUs, its implementation is more complicated than spatial blocking [28]. For example, for every GPU, it is necessary to calculate and ensure that a thread block has enough shared and local memory to save data for calculating a group of time/computation steps. Consequently, examples of stencil computation in GPUs and temporal blocking are rare. To apply this technique to CA (especially to the most complex ones), automatic code generators should be used. The two environments that do this in GPUs efficiently are STENCILGEN [29] and AN5D [30]. In addition, the AN5D currently has the best performance results in stencil computing and GPUs and it is publicly available for use. AN5D improves the 3.5D blocking algorithm obtaining the best performance in stencil computation for both single and double precision floating point. Only a 𝑝𝑟𝑎𝑔𝑚𝑎 directive is necessary to indicate which loop must be translated into Cuda code. In AN5D, the C code that is translated into Cuda code must be composed only of three ‘for’ loops and a variable assignment. It is not possible to include selective structures (if-else) or previous instructions that have dependencies on the cell values. AN5D must Simulation Modelling Practice and Theory 118 (2022) 102519 8 D. Cagigas-Muñiz et al. Fig. 6. Basic example of overlapped tiling (temporal blocking) in a GPU. resolve dependencies at compiling time. For all these reasons, look-up tables (see Section 5.1) must be used too. Therefore, the C sequential code must be codified as shown below, so that AN5D can be compiled to cuda: ... #pragma scop for ( int t = 0; t < TIMESTEPS ; t++) for ( int i = 1; i < SIZE −1; i ++) for ( int j = 1; j < SIZE −1; j ++) { grid [( t +1)%2][ i ][ j ] = lookup_table [ grid [ t %2][ i ][ j ]] [ ( grid [ t%2][i −1][ j ] + grid [ t %2][ i +1][ j ] + Simulation Modelling Practice and Theory 118 (2022) 102519 9 D. Cagigas-Muñiz et al. Table 1 Mean execution time and standard deviation (left and right values between parenthesis in each table cell) of 15 independent runs in milliseconds of the Game of Life (GOL) cellular automaton for different GPU implementations: baseline, baseline using a look-up table to code rules, temporal blocking using AN5D framework (with default options), packet coding using 32 bits per cell (4 subcells of 8 bits) and packet coding using 64 bits per cell (8 subcells of 8 bits). Best case for each algorithm/technique and each grid size is presented in bold. Grid size Baseline Cuda implementation Look-up table AN5D (Temporal blocking) 32 bits packet coding (4 subcells per cell) 64 bits packet coding (8 subcells per cell) 256 6.2 (0.0) 5.8 (0.0) 25.7 (0.1) 3.7 (0.0) 3.9 (0.0) 512 20.9 (0.4) 20.4 (0.3) 34.9 (0.2) 8.9 (0.0) 10.5 (0.2) 1024 88.5 (9.7) 87.4 (9.7) 82.3 (7.5) 39.6 (0.3) 38.6 (0.8) 2048 325.2 (11.6) 305.0 (0.9) 264.8 (1.8) 128.1 (8.0) 114.9 (9.3) 4096 1299.2 (2.6) 1234.4 (1.2) 1036.0 (8.8) 488.5 (1.5) 436.0 (1.4) 8192 5300.5 (16.6) 5052.2 (16.4) 4228.5 (18.2) 1977.6 (8.6) 1770.2 (10.8) 16384 21303.0 (60.8) 20324.4 (58.6) 17066.8 (71.8) 7960.2 (31.4) 6996.2 (27.2) grid [ t %2][ i ][ j −1] + grid [ t %2][ i ][ j +1] + grid [ t %2][i −1][j −1] + grid [ t %2][i −1][ j +1] + grid [ t %2][ i +1][j −1] + grid [ t %2][ i +1][ j +1]) ] ; } #pragma endscop ... The variable assignment inside the three loops calculates the new cell state based on neighbour cells. The main data structure used should be only a 3 dimensional array. First array dimension contains the current and new grids. The grids swap their roles in each time step. Second and third array dimensions contain the rows and columns. 6. Experimental results In this section, the algorithms and techniques proposed before are verified experimentally. Four representative CA that have been proposed in the scientific bibliography are selected. Two CA have Moore neighbourhood and two have Von Neuman neighbourhood. The number of states ranges from 2 to 15. Rules are diverse and the NVIDIA profiler indicates memory bandwidth overhead when using the standard baseline implementation version. In the next subsections, the results obtained for each cellular automaton are described. For all experiments in this paper, a NVIDIA GTX 1650 TI graphics card was used in a computer with an AMD Ryzen 5 3600 processor. Linux-Mint 20.04 operating system, NVCC 10.1, and GCC 9.3 compilers were part of the software configuration. Cuda programming language was used to code CA. Every test of this article was executed 15 times and the average and standard deviation are presented. Experiments (Cuda implementations) use different CA grid sizes. From 256 ×256 cells to 16384 ×16384 cells. Results are detailed in tables instead of graphics to appreciate the performance differences in small grid sizes. The baseline Cuda implementation result is always placed at the first table column. This is the reference implementation when comparing with the other results. Only the CA main Cuda kernel is measured using the NVIDIA Profiler in each case. Initial procedures like setting up the original grid states, are not taken into account. The main CA kernel computes the time steps. It represents between 72% (256 ×256 grid) and 99% (16384 ×16384 grid) of the total execution time for the baseline Cuda implementation. For all experiments, 1024 time steps were used. The source code of all experiments is available at [31]. 6.1. Game of life The Game of Life (GoL) cellular automaton was the first proposed cellular automaton. This simple cellular automaton is ideal to make a first test on the algorithms and techniques proposed in the sections before. In Table 1 there is a summary of the results obtained. Results in Table 1 include GoL Cuda implementation using only LUTs, AN5D framework, packet coding of 32 bits per cell, and packet coding of 64 bits per cell. The average execution times in the baseline and LUTs versions are very similar. The AN5D version has better execution time but only for 1024 ×1024 grid size and above. In the case of 16384 ×16384 GoL grid (the biggest grid), there is a speed up of 25% when comparing with GoL baseline execution time. Packet coding implementations present the best results: a speed-up above 300% is achieved for the largest CA. However, there is almost a 14% time improvement when using the 64 bit packet coding algorithm (8 subcells per cell) instead of the 32 bit version (4 subcells per cell). Therefore, the next CA will use this packet coding version. A Kruskal–Wallis test [32] was performed to determine whether or not there was a statistically significant difference between the median run times of the different techniques or algorithms used. Each technique or algorithm was applied to 7 different grid sizes (from 256 to 16384) as shown in Table 1. Simulation Modelling Practice and Theory 118 (2022) 102519 16 D. Cagigas-Muñiz et al. References [1] S.F. Judice, Lattice gas cellular automata for fluid simulation, in: Encyclopedia of Computer Graphics and Games, Springer International Publishing, Cham, 2018, pp. 1–8, http://dx.doi.org/10.1007/978-3-319-08234-9_184-1. [2] B. Arca, T. Ghisu, G.A. Trunfio, GPU-accelerated multi-objective optimization of fuel treatments for mitigating wildfire hazard, J. Comput. Sci. 11 (2015) 258–268, http://dx.doi.org/10.1016/j.jocs.2015.08.009. [3] R. Lubas, J. Was, J. Porzycki, Cellular automata as the basis of effective and realistic agent-based models of crowd behavior, J. Supercomput. 72 (6) (2016) 2170–2196, http://dx.doi.org/10.1007/s11227-016-1718-7. [4] J. Kroc, F. Jiménez-Morales, J.L. Guisado, M.C. Lemos, J. Tkáč, Building efficient computational cellular automata models of complex systems: background, applications, results, software, and pathologies, Adv. Complex Syst. 22 (05) (2019) 1950013. [5] K.R. Tubbs, F.T.-C. Tsai, GPU accelerated lattice Boltzmann model for shallow water flow and mass transport, Internat. J. Numer. Methods Engrg. 86 (3) (2011) 316–334, http://dx.doi.org/10.1002/nme.3066, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.3066. [6] A.G. Salguero, A.J. Tomeu-Hardasmal, M.I. Capel, Dynamic load balancing strategy for parallel tumor growth simulations, J. Integr. Bioinform. 16 (1) (2019) 20180066, http://dx.doi.org/10.1515/jib-2018-0066. [7] M. Sitko, K. Banas, L. Madej, Scaling scientific cellular automata microstructure evolution model of static recrystallization toward practical industrial calculations, Materials 14 (2021) (2021) http://dx.doi.org/10.3390/ma14154082. [8] B. Jelinek, M. Eshraghi, S. Felicelli, J.F. Peters, Large-scale parallel lattice Boltzmann-cellular automaton model of two-dimensional dendritic growth, Comput. Phys. Comm. 185 (3) (2014) 939–947, http://dx.doi.org/10.1016/j.cpc.2013.09.013. [9] C. Xia, H. Wang, A. Zhang, W. Zhang, A high-performance cellular automata model for urban simulation based on vectorization and parallel computing technology, Int. J. Geogr. Inf. Sci. 32 (2) (2018) 399–424, http://dx.doi.org/10.1080/13658816.2017.1390118. [10] A. Kerr, G. Diamos, S. Yalamanchili, A characterization and analysis of PTX kernels, in: 2009 IEEE International Symposium on Workload Characterization, IISWC, 2009, pp. 3–12, http://dx.doi.org/10.1109/IISWC.2009.5306801. [11] M.J. Gibson, E.C. Keedwell, D.A. Savić, An investigation of the efficient implementation of cellular automata on multi-core CPU and GPU hardware, J. Parallel Distrib. Comput. 77 (2015) 11–25, http://dx.doi.org/10.1016/j.jpdc.2014.10.011. [12] S. Rybacki, J. Himmelspach, A. Uhrmacher, CPU and GPU based simulation of cellular automata - a performance comparison, in: Proceedings of the 1st SIMUL, 2009, pp. 62–67. [13] E. Millán, P. Martínez, G. Gil Costa, M. Piccoli, A. Printista, C. Bederian, C. García Garino, E. Bringa, in: A. De Giusti (Ed.), Parallel implementation of a cellular automata in a hybrid CPU/GPU environment, XVIII Congreso Argentino de Ciencias de la Computación, 2013, pp. 184–193. [14] E.R. Berlekamp, J.H. Conway, R.K. Guy, Winning Ways for Your Mathematical Plays, second ed., A K Peters/CRC Press, New York, USA, 2001. [15] C. Daniel, F. Diaz-del Rio, M. López-Torres, F. Jiménez-Morales, J.L. Guisado, Developing efficient discrete simulations on multicore and GPU architectures, Electronics 9 (2020) 189, http://dx.doi.org/10.3390/electronics9010189. [16] E. Nicolas, N. Wolovick, F. Piccoli, C. Garcia Garino, E. Bringa, Performance analysis and comparison of cellular automata GPU implementations, Cluster Comput. 20 (2017) (2017) http://dx.doi.org/10.1007/s10586-017-0850-3. [17] W.-m.W. Hwu, GPU Computing Gems Jade Edition, first ed., Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 2011. [18] F. Diaz-del Rio, D. Cagigas-Muñiz, J.-L. Guisado-Lizar, J.-L. Sevillano-Ramos, Efficient parallel implementation of cellular automata and stencil computations in current processors, in: P. Nicopolitidis, S. Mistra, L. Yang, B. Zeigler, Z. Ning (Eds.), Advances in Computing, Informatics, Networking and Cybersecurity - a Book Honoring Prof. Mohammad S. Obaidat’s Significant Scientific Contributions, Springer-Nature, 2022, pp. 1–29. [19] G. Ofenbeck, R. Steinmann, V. Caparros, D.G. Spampinato, M. Püschel, Applying the roofline model, in: 2014 IEEE International Symposium on Performance Analysis of Systems and Software, ISPASS, 2014, pp. 76–85, http://dx.doi.org/10.1109/ISPASS.2014.6844463. [20] A. Simpson, Oak ridge leadership computing facility, URL https://github.com/olcf/game_of_life_tutorials/tree/master/CUDA. [21] G. Oxman, S. Weiss, Y. Be’ery, Computational methods for Conway’s game of life cellular automaton, J. Comput. Sci. 5 (2013) (2013) http://dx.doi.org/ 10.1016/j.jocs.2013.07.005. [22] D.B. Kirk, W.-m.W. Hwu, Programming Massively Parallel Processors: A Hands-on Approach, Morgan Kaufmann Publishers, Burlington, MA, 2010. [23] P.S. Rawat, M. Vaidya, A. Sukumaran-Rajam, A. Rountev, L. Pouchet, P. Sadayappan, On optimizing complex stencils on GPUs, in: 2019 IEEE International Parallel and Distributed Processing Symposium, IPDPS, 2019, pp. 641–652. [24] A. Schäfer, D. Fey, High performance stencil code algorithms for GPGPUs, Procedia Comput. Sci. 4 (2011) 2027–2036, http://dx.doi.org/10.1016/j.procs. 2011.04.221. [25] J. Holewinski, L.-N. Pouchet, P. Sadayappan, High-performance code generation for stencil computations on GPU architectures, in: Proceedings of the 26th ACM International Conference on Supercomputing, New York, NY, USA, 2012, pp. 311–320, http://dx.doi.org/10.1145/2304576.2304619. [26] P. Rawat, Optimization of Stencil Computations on GPUs (Electronic Thesis Or Dissertation), Ohio State University, 2018. [27] A.D. Nguyen, N. Satish, J. Chhugani, C. Kim, P. Dubey, 3.5-D blocking optimization for stencil computations on modern CPUs and GPUs, in: SC, IEEE, 2010, pp. 1–13. [28] K. Hou, H. Wang, W.-c. Feng, Gpu-UniCache: Automatic code generation of spatial blocking for stencils on GPUs, in: Proceedings of the Computing Frontiers Conference, in: CF’17, Association for Computing Machinery, New York, NY, USA, 2017, pp. 107–116, http://dx.doi.org/10.1145/3075564.3075583. [29] P.S. Rawat, M. Vaidya, A. Sukumaran-Rajam, M. Ravishankar, V. Grover, A. Rountev, L.-N. Pouchet, P. Sadayappan, Domain-specific optimization and generation of high-performance GPU code for stencil computations, Proc. IEEE 106 (11) (2018) 1902–1920. [30] K. Matsumura, H. Zohouri, M. Wahib, T. Endo, S. Matsuoka, AN5D: automated stencil framework for high-degree temporal blocking on GPUs, in: International Symposium on Code Generation and Optimization, 2020, pp. 199–211, http://dx.doi.org/10.1145/3368826.3377904. [31] D.C.-M. niz, Cellular automata software repository, URL https://github.com/dcagigas/GPU-Cellular-Automata. [32] E. Ostertagova, O. Ostertag, J. Kováč, Methodology and application of the Kruskal-Wallis test, Appl. Mech. Mater. 611 (2014) 115–120. [33] D. Rey, M. Neuhäuser, Wilcoxon-signed-rank test, international encyclopedia of statistical science, in: International Encyclopedia of Statistical Science, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011, pp. 1658–1659, http://dx.doi.org/10.1007/978-3-642-04898-2_616. [34] J.G. Freire, C.C. DaCamara, Using cellular automata to simulate wildfire propagation and to assist in fire management, Nat. Hazards Earth Syst. Sci. 19 (1) (2019) 169–179, http://dx.doi.org/10.5194/nhess-19-169-2019, URL https://nhess.copernicus.org/articles/19/169/2019/. [35] Y. Zhao, D. Geng, Simulation of forest fire occurrence and spread based on cellular automata model, in: ICAIIS 2021: 2021 2nd International Conference on Artificial Intelligence and Information Systems, Chongqing, China, May 28 - 30, 2021, ACM, 2021, pp. 304:1–304:6, http://dx.doi.org/10.1145/3469213. 3471332. [36] L. Hugo, P. Hugo, P. Thomas, AutoCelle in C++, URL https://github.com/hugofloter/CelularAutomaton. [37] D. Griffeath, Self-organizing two-dimensional cellular automata: 10 still frames, in: Designing Beauty: The Art of Cellular Automata, Springer International Publishing, Cham, 2016, pp. 1–12, http://dx.doi.org/10.1007/978-3-319-27270-2_1. [38] K. Kwak, Y. Baryshnikov, E. Coffman, Cyclic cellular automata: A tool for self-organizing sleep scheduling in sensor networks, in: Proceedings - 2008 International Conference on Information Processing in Sensor Networks, IPSN 2008, 2008, pp. 535–536, http://dx.doi.org/10.1109/IPSN.2008.69. [39] R. González-García, G. Castanon, H.E. Hernández Figueroa, 2D photonic crystal complete band gap search using a cyclic cellular automaton refination, Photon. Nanostruct.: Fundam. Appl. 12 (2014) (2014) http://dx.doi.org/10.1016/j.photonics.2014.09.003. [40] V. Gladkikh, A. Nigay, Wireworld++: A cellular automaton for simulation of nonplanar digital electronic circuits, Complex Systems 27 (2018) (2018). [41] C. Luo, R. Suda, A performance and energy consumption analytical model for GPU, in: 2011 IEEE Ninth International Conference on Dependable, Autonomic and Secure Computing, 2011, pp. 658–665, http://dx.doi.org/10.1109/DASC.2011.117.