Simulating the behavior of the human brain on GPUS

Simulating the behavior of the Human Brain on GPUs Pedro Valero-Lara 1,* , Ivan Martı ´nez-Pe ´rez 1 , Rau ¨l Sirvent 1 , Antonio J. Pen ˜a 1 , Xavier Martorell 2 , and Jesu ´s Labarta 2 1 Barcelona Supercomputing Center (BSC), C/Jordi Girona, 29, 08034 Barcelona, Spain 2 Universidad Polite `cnica de Catalun ˜a, Carrer de Jordi Girona, 1, 3, 08034 Barcelona, Spain Received: 23 February 2018 / Accepted: 11 September 2018 Abstract. The simulation of the behavior of the Human Brain is one of the most important challenges in computing today. The main problem consists of finding efficient ways to manipulate and compute the huge volume of data that this kind of simulations need, using the current technology. In this sense, this work is focused on one of the main steps of such simulation, which consists of computing the Voltage on neurons’ morphology. This is carried out using the Hines Algorithm and, although this algorithm is the optimum method in terms of number of operations, it is in need of non-trivial modifications to be efficiently parallelized on GPUs. We proposed several optimizations to accelerate this algorithm on GPU-based architectures, exploring the limitations of both, method and architecture, to be able to solve efficiently a high number of Hines systems (neurons). Each of the optimizations are deeply analyzed and described. Two different approaches are studied, one for monomorphology simulations (batch of neurons with the same shape) and one for multi-morphology simulations (batch of neurons where every neuron has a different shape). In mono-morphology simulations we obtain a good performance using just a single kernel to compute all the neurons. However this turns out to be inefficient on multi-morphology simulations. Unlike the previous scenario, in multi-morphology simulations a much more complex implementation is necessary to obtain a good performance. In this case, we must execute more than one single GPU kernel. In every execution (kernel call) one specific part of the batch of the neurons is solved. These parts can be seen as multiple and independent tridiagonal systems. Although the present paper is focused on the simulation of the behavior of the Human Brain, some of these techniques, in particular those related to the solving of tridiagonal systems, can be also used for multiple oil and gas simulations. Our studies have proven that the optimizations proposed in the present work can achieve high performance on those computations with a high number of neurons, being our GPU implementations about 4·and 8·faster than the OpenMP multicore implementation (16 cores), using one and two NVIDIA K80 GPUs respectively. Also, it is important to highlight that these optimizations can continue scaling, even when dealing with a very high number of neurons. 1 Motivation Today, we can find multiple initiatives that attempt to simulate the behavior of the Human Brain by computer [1–3]. This is one of the most important challenges in the recent history of computing with a large number of practical applications. The main constraint is being able to simulate efficiently a huge number of neurons using the current computer technology. One of the most efficient ways in which the scientific community attempts to simulate the behavior of the Human Brain consists of computing the next three major steps [4]: The computing of (1) the Voltage on neuron morphology, (2) the synaptic elements in each of the neurons and (3) the connectivity between the neurons. In this work, we focus on the first step which is one of the most time consuming steps of the simulation. Also, it is strongly linked with the rest of steps. All these steps must be carried out on each of the neurons. The Human Brain is composed by about 11 billion of neurons, which are completely different among them in size and shape. The standard algorithm used to compute the Voltage on neurons’ morphology is the Hines algorithm [5], which is based on the Thomas algorithm [6], that solves tridiagonal systems. Although the use of GPUs to compute the Thomas algorithm has been deeply studied [7–11], the differences among these two algorithms, Hines and Thomas, make us impossible to use the last one, as this cannot deal with the sparsity of the Hines matrix. The solving of one Hines system can be also seen as a set of independent and non-independent triangular systems, which could be solved by using the Thomas algorithm. Previous works [12] have explored the use of other algorithms based on the Stone’s method [13]. Unlike Thomas algorithm, this method is parallel. However, it is in need * Corresponding author: [email protected] This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018) Available online at: P. Valero-Lara et al., published by IFP Energies nouvelles, 2018 www.ogst.ifpenergiesnouvelles.fr https://doi.org/10.2516/ogst/2018061 REGULAR ARTICLEREGULAR ARTICLE Numerical methods and HPC A. Anciaux-Sedrakian and Q.H. Tran (Guest editors) of a higher number of operations (20nlog 2n)withrespect to the (8n) operations of the Thomas algorithm to solve one single system of size n. Also, the use of parallel methods presents some additional drawbacks to be dealt with. For instance, it would be difficult to compute those neurons that compromise a size bigger than the maximum number of threads per CUDA block (1024) or shared memory (48 KB). Unlike the work presented in [12], where a relatively low number of neurons (128) is computed using single precision operations, in this work we are able to execute a very high number of neurons (up to hundreds of thousands) using double precision operations. We have used the Hines algorithm, which is the optimum method in terms of number of operations, avoiding high expensive computational operations, such as synchronizations and atomic accesses. Our code is able to compute a high number of systems (neurons) of any size in one call (CUDA kernel), using one thread per Hines system instead of one CUDA block per system. Although multiple works have explored the use of GPUs to compute multiple independent problems in parallel without transforming the data layout [14–17], the particular characteristics of the sparsity of the Hines matrices force us to modify the data layout to efficiently exploit the memory hierarchy of the GPUs (coalescing accesses to GPU memory). The present work extends the previously published work [18] with additional contributions. This work includes a complete new approach to deal with one of the most important challenges in the simulation of the Human Brain, that is, dealing with simulations which involve neurons with different morphologies (multi-morphology simulations). To deal with this particular scenario, we must compute parts of the batch of neurons separately. These parts can be seen as multiple and independent tridiagonal systems. While the present paper is focused on the simulation of the behavior of the Human Brain, some of these techniques, in particular those related to the solving of tridiagonal systems, can be also used for multiple oil and gas simulations. 1 This article is structured as follows: Section 2 briefly introduces the physical problem at hand and the general numerical framework that has been selected to cope with it: Hines algorithm. In Section 3 we present the specific parallel features for the resolution of multiple Hines systems for mono-morphology simulations, as well as the parallel strategies envisaged to optimally enhance the performance. Section 4 shows the strategies proposed for dealing with the challenges presented in the multi-morphology simulations. The state-of-the-art references and the differences between these and the present work are presented in Section 5. Finally, the conclusions are outlined in Section 6. 2 Hines algorithm In this section, we describe the numerical framework behind the computation of the Voltage on neurons morphology. It follows the next general form: CoV otþI¼fo oxgoV ox  ð1Þ where fand gare functions on x-dimension and the current Iand capacitance C[4] depend on the voltage V. Discretizing the previous equation on a given morphology we obtain a system that has to be solved every time-step. This system must be solved at each point: aiVnþ1 iþ1þdiVnþ1 iþbiVnþ1 i1¼rið2Þ where the coefficients of the matrix are defined as follows: Upperdiagonal : ai¼figiþ1 2 2D2 x Lowerdiagonal : bi¼figiþ1 2 2D2 x Diagonal : di¼Ci DtðaiþbiÞ rhs :ri¼Ci DtVn iIaiðVn i1Vn iÞbiðVn iþ1Vn iÞ a i and b i are constant in the time, and they are computed once at start up. Otherwise, the diagonal (d) and righthand-side (rhs) coefficients are updated every time-step when solving the system. The discretization above explained is extended to include branching, where the spatial domain (neuron morphology) is composed of a series of one-dimension sections that are joined at branch points according to the neuron morphology. For the sake of clarity, we illustrate a simple example of a neuron morphology in Figure 1. It is important to note that the graph formed by the neuron morphology is an acyclic graph, i.e. it has no loops. The nodes are numbered using a scheme that gives the matrix sparsity structure that allows to solve the system in linear time. To describe the sparsity of the matrix from the numbering used, we need an array (p i i2[2:n]) which stores the parent indexes of each node. The pattern of the matrix which illustrates the morphology shown above is graphically illustrated in Figure 1. The Hines matrices feature the following properties: they are symmetric, the diagonal coefficients are all nonzero and per each off-diagonal element, there is one off-diagonal element in the corresponding row and column (see row/column 7, 12, 17 and 22 in Fig. 1). Given the aforementioned properties, the Hines systems (Ax =b) can be efficiently solved by using an algorithm similar to Thomas algorithm for solving tri-diagonal systems. This algorithm, called Hines algorithm, is almost identical to the Thomas algorithm except by the sparsity pattern given by the morphology of the neurons whose pattern is stored by the pvector. An example of the sequential code used to implement the Hines algorithm is illustrated in pseudo-code in Algorithm 1. 1 https://tridiagonal.com/oil-and-gas P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018)2 3 Implementation of batch mono-morphology Hines on GPUs Efficient memory management is critical to achieve a good performance, but even much more on those architectures based on high throughput and high memory latency, such as GPUs. In this sense, we focus on presenting the different data layouts proposed and analyze the impact of these on the overall performance. Three different data layouts were explored: Flat, Full-Interleaved and Block-Interleaved. While the Flat data layout consists of storing all the elements of each of the systems in contiguous memory locations, inthe Full-Interleaved data layout, we start by storing the first elements of each of the systems in contiguous memory locations, after that we store the set of the second elements, and so on until the last element. Similarly to the Full-Interleaved data layout, the Block-Interleaved data layout divides the set of systems into groups of systems of a given size (BS), whose elements are stored in memory by using the strategy followed by the Full-Interleaved approach. For the sake of clarity, Figure 2 illustrates a simple example composed by four different Hines systems of three elements each. Please, note that we only illustrate one vector per system in Figure 2, but in the real scenario we would have four vectors per Hines system (Pseudocode 1) on which the strategies above described are carried out. As widely known, one of the most important requirements to achieve a good performance on NVIDIA GPUs is to have contiguous threads accessing contiguous memory locations Fig. 1. Example of a neuron morphology and its numbering (left-top and bottom) and sparsity pattern corresponding to the numbering followed (top-right) [19]. Size of the system Set of vectors Set of first elements Set of second elements Set of third elements Second blockFirst block Fig. 2. Example of the different data layouts proposed, Flat (top), Full-Interleaved (center), Block-Interleaved (bottom) with aBS equal to 2, for four Hines systems of three elements each. Dotted lines represent the jumps in memory carried out by the first thread/system. Algorithm 1 Hines algorithm. 1 void solveHines(double*u, double *l, double *d, 2 double *rhs, int *p, int cellS ize) 3//u^upper vector, l ^lower vector 4 inti; 5 double factor; 6 // Backward Sweep 7 for i =cellS ize - 1 ^0do 8 factor =u[i] /d[i]; 9 d[p[i]] -=factor xl[i]; 10 rhs[p[i]] -=factor xrhs[i]; 11 end for 12 [0 ] d [ / = [0 ] h sr 13 // Forward Sweep 14 for i =1^cellS ize - 1 do 15 rhs[i] -=l[i] xrhs[p[i]]; 16 rhs[i] /= d[i]; 17 end for P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018) 3 (coalescing memory accesses). This is the main motivation behind the proposal of the different data layouts. Our GPU implementation consists of using one thread perHinessystem.AscommentedinSection 1, we decided to explore this approach to avoid dealing with atomic accesses and synchronizations, as well as to be able to execute a very high number of neurons of any size. Using the Flat data layout we cannot exploit coalescence; however by interleaving (Full-Interleaved data layout) the elements of the vectors (u,l,d,rhs and pin Pseudocode 1), contiguous threads access to contiguous memory locations. Although we exploit coalescence in memory accesses by using this approach, the threads have to jump in memory as many elements as the number of systems to access the next element of the vector(s) (dotted lines in Fig. 2). This could cause an inefficient use of the memory hierarchy. This is why we study an additional approach, the called BlockInterleaved data layout. Using this approach we reduce the number of elements among consecutive elements of the same system, and hence the jumps in memory are not as big as in the previous approach (Full-Interleaved), while keeping the coalesced memory accesses. Also, the use of the Block-Interleaved data layout can take better advantage of the growing importance of the bigger and bigger cache memories in the memory hierarchy of the current and upcoming GPU architectures. 3.1 Implementation based on Shared Memory Unlike the previous approaches, here we explore the use of shared memory for our target application. The shared memory is much faster than the global memory, but it presents some important constraints to deal with. This memory is useful when the same data can be reused either by the same thread or by other thread of the same block of threads (CUDA block). Also, it is small (up to 48 KB in the architecture used) and its use hinders the exchange among blocks of threads by the CUDA scheduler to overlap accesses to global memory with computation. As we can see in Pseudocode 1, in our problem the elements of the vectors a, d, b, rhs and pare reused in the Forward Sweep after computing the Backward Sweep. However, the granularity used (1 thread per system) and the limit of the shared memory (48 KB) prevent from storing all the vectors in shared memory. To be able to use shared memory we have to use the Block-Interleaved data layout. The number of systems to be grouped (BS)isimposedby the size of the shared memory. In order to address the limitation of shared memory, we only store the rhs vector, as this is the vector on which more accesses are carried out. In this sense, the more systems are packed in shared memory, the more accesses to shared memory are carried out. 3.2 Performance analysis For the experiments, we have used a heterogeneous node 2 composed of 2·Intel Xeon E5-2630v3 (Haswell) with 8 cores and 20 MB L3 cache each, and 2·K80 NVIDIA GPU (Kepler) with a total of 4992 cores and 24 GB GDDR5 of global memory each. Each K80 is composed of 2·logic GPUs similar to K40. This node is a Linux (Red Hat 4.4.716) machine, on which we have used the next configuration (compilers version and flags): gcc 4.4.7, nvcc (CUDA) 7.5, -O3, -fopenmp, -arch=sm_37. The code evaluated in this section is available in a public access repository. 3 To evaluate the different implementations described in the previous section, we have used real configurations (neurons’ morphologies). 4 In particular, six different neurons were used, which can be divided into six different categories regarding their sizes and number of branches. More details are described in Table 1. We have considered these six different morphologies, as a wide range of the neurons fall into the chosen morphologies. In this Section, five different implementations are analyzed. One is based on OpenMP Multicore using the Flat data layout (Sect. 3), which makes use of an OpenMP pragma (#pragma omp for) on the top of the for loop which goes over the different independent Hines systems to distribute blocks of systems over the available cores. The rest of implementations are based on GPU. Basically, we have one implementation per each of the data layout described: Flat, Full-Interleaved and Block-Interleaved. Additionally, we study the use of shared memory (Block-Shared) over Block-Interleaved. There are multiple different configurations regarding the block-size (BS) and CUDA block for the last two scenarios (Block-Interleaved and BlockShared). For sake of clarity we focus on one of the possible Table 1. Summary of the neurons used. Name Size #Branches Code name Neuron ID small-low 76 7 299-DG-IN-Neuron2 NMO_00076 small-high 76 29 202-2-19nj NMO_00076 medium-low 305 30 59D-40X NMO_00302 medium-high 319 157 Culture-9-5 NMO_00319 big-low 695 66 28-2-2 NMO_00695 big-high 691 341 HSE-fluoro02 NMO_00691 2 MinoTauro, https://www.bsc.es/ca/innovation-and-services/ supercomputers-and-facilities/minotauro 3 BSC-GitLab, https://pm.bsc.es/gitlab/imartin1/cuHinesBatch 4 http://www.neuromorpho.org/ P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018)4 test cases to evaluate these two approaches. The benefit shown for these two implementations is similar to the rest of test-cases. First, we evaluate the Multicore,Flat and FullInterleaved for 256, 2560, 25 600 and 256 000 mediumhigh (Tab. 1)neurons(Fig. 3). On those test cases that do not compromise a high number of neurons (256 and 2560), Multicore obtains better performance than the GPU-based implementations. This is mainly because of the parallelism of these tests, which is not enough to saturate GPU and this cannot reduce the impact of the high latency by overlapping execution and memory accesses. The use of multicore (16 cores and 2 sockets) supposes a speedup (over sequential execution) about 2 for 256 neurons and about 6 for 256 000 neurons. As shown, Flat is not able to scale, even on those test-cases that involve a high number of neurons, being even slower than multicore execution, achieving a maximum speedup of about 2. This is because of the memory access pattern which cannot exploit coalescing (contiguous threads access to contiguous memory locations). On the other hand, Full-Interleaved turns out as the best choice, being faster than Multicore and Flat, when dealing with a high number of neurons (25 600 and 256 000). Unlike Flat, Full-Interleaved takes advantage of coalescing when accessing to global memory. As expected, this has an impressive impact on performance, being Full-Interleaved about 20·and 25·faster than sequential code when computing 25 600 and 256 000 neurons respectively on one K80 GPU. The use of multiple GPUs is only beneficial on those test cases with an enough computational load where a high number of neurons must be computed (25 600 and 256 000 neurons), with an extra benefit close to the ideal scaling (about 1.9·faster than using one K80 GPU). As it is not possible to have control on the CUDA scheduler, we have explored a high number of different combinations regarding block-size (BS) for the Block-Interleaved approach (Sect. 3). For the sake of clarity, and given the huge number of different possible test-cases, we have focused on one particular scenario. It consists of computing 256 000 medium-high neurons using different block sizes (BS) and fixing the size of the CUDA block (number of threadsperblock).Thisisacharacteristiccaseamongthe tests carried out, as the features (sizes and number of branches) of the morphology used is in between of the other two morphologies. As shown in Figure 4a, some of the cases are slightly better than the Full-Interleaved approach, being about a 2% faster. Next we analyze the performance of the Block-Shared implementation. We focus on the same scenario used for the Block-Interleaved.Figure 4b graphically illustrates the performance achieved by the Block-Shared and the other approaches. Although using shared memory is better than the performance achieved by the Flat approach, it is much smaller than the Full-Interleaved counterpart. For this particular scenario (medium-high morphology), a very low number of systems saturate the capacity of the shared memory (48 KB). Also, the data reuse is low using onethread per Hines system. These drawbacks do not allow to achieve a better performance when the shared memory is used. Finally, we evaluate the impact on performance of the particularities of each of the morphologies (Tab. 1). The performance achieved by the Flat is not included as it was proven to be very inefficient. As shown in Figure 5, both approaches, Multicore and Full-Interleaved, show a similar trend in performance independently of the neurons’ morphology. In particular, the peak speedup achieved on the different morphologies does not vary significantly (47·–55·). After comparing the performance achieved by Multicore and GPU, now we focus on evaluating the efficiency of our GPU implementation. To do that, we make use of nvprof. 5 We do not obtain very different results depending on the input (number and shape of neurons). In all cases, we obtain more than 99% efficiency (sm_efficiency), as well as a bandwidth (Global Load Throughput) close to 160 GB/s, being the theoretical peak equal to 240 GB/s and the effective about the bandwidth achieved by our implementation. As most of the GPU applications, our implementation is memory bound and this is reflected by a low occupancy (about 24%). 3.3 Remarks Block-Interleaved is positioned as the fastest approach against the others when dealing with a high number of neurons. However this implementation is difficult to tune. It is not possible to know the best configuration in advance. In contrast, Full-Interleaved is almost as fast as BlockInterleaved and is not in need of being tuned a priori. It is important to highlight that both approaches require to modify the data layout by interleaving the elements of the vectors. This preprocessing compromises an irrelevant cost with respect to the whole process, as for our target application (the simulation of the Human Brain), this is 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines Systems Multicore Flat(one logic GPU) Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) Fig. 3. Performance (speedup over sequential execution) achieved by Multicore (16 cores, 2 sockets) and the GPU-based approaches, Flat and Full-Interleaved (using different number of GPUs), using medium-high neurons. 5 nvprof -m achieved_occupancy,sm_efficiency,gld_throughput, gst_throughput,gld_efficiency,gst_efficiency ./run P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018) 5 carried out just once at the very beginning of the simulation. Using Full-Interleaved we obtain a similar behavior in terms of performance when different morphologies are considered. This is particularly interesting to evaluate the scalability and robustness of the implementation. As overview, while Multicore (16 cores and 2 sockets) gives us a maximum speedup against sequential execution of about 6·, reporting a peak speedup of about 55·,using2·K80 NVIDIA GPUs. 4 Implementation of batch multi-morphology Hines on GPUS In this section, we analyze the performance for the simulation of the behavior of the Human Brain using different morphologies (multi-morphology) in parallel. This study is the most important contribution of the present work. Although the simulations that involve mono-morphology (the same neuron replicated) scenarios can also be useful, in real world cases, the neurons are completely different among them. First, we evaluate the optimizations above described for this particular case. Figure 6 graphically illustrates the performance achieved for multi-morphology simulations using the strategies described in the previous section for the mono-morphology approach.Todothis, we have used two different test cases. In both cases we use the same size to evaluate only the influence on performance for computing multi-morphology simulations. We generate different Hines matrices with a different ratio (% of branches with respect to the size of branches). For the sake of comparison, we also include those cases for mono-morphology simulations using the same morphology, size and branches ratio for the batch of neurons (Mono in Fig. 6). As shown (Fig. 6), when dealing with multimorphology cases, we found an important fall in performance with respect to mono-morphology simulations. The most important cause of this behavior is given by the lack of coalescing memory accesses. When computing different morphologies in parallel, the off-diagonal elements (see Fig. 1 and Pseudocode 1) of these Hines matrices (neurons) are located in different positions from one to other, which makes difficult that consecutive CUDA threads of the same CUDA block access to contiguous memory positions causing an important fall in performance. To minimize this problem, we propose a different approach which consists of computing each of the branch levels separately. Each branch can be seen as a tridiagonal system, so those branches of the same level could be computed simultaneously in one CUDA kernel. For the sake of clarity, Figure 7 shows a simple diagram with this idea. In the rest of this section, we focus on the implementation of a kernel, which makes use of some of the ideas previously presented, but to solve tridiagonal systems instead of Hines systems. At the end of this section, we evaluate the impact of this idea for multi-morphology simulations. 4.1 Tridiagonal linear systems The state-of-the-art method to solve tridiagonal systems is the called Thomas algorithm [8], which a specialized application of the Gaussian elimination that takes into account the tridiagonal structure of the system. It consists of two stages, commonly denoted as forward elimination and backward substitution. Given a linear Au =ysystem, where Ais a tridiagonal matrix: A¼ b1c10 a2b2c2 :: : :: : an1bn1cn1 anbn 2 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 5 : 12.4 12.6 12.8 13 13.2 Full-Inter 32 64 128 256 512 SpeedUp Block(Interleaved)-Size 2 4 6 8 10 12 14 Multi Flat Full-Inter Shared SpeedUp Approaches (a) (b) Fig. 4. (a) Performance (speedup over sequential execution) achieved by the Block-Interleaved approach for multiple BS (32, 64, 128, 256, 512) for a CUDA Block size equal to 128. (b) Performance (speedup over sequential execution) achieved by the Block-Shared implementation, Flat, Full-Interleaved (Full-Inter) and Multicore (Multi) using 16 cores. The test-case consisted of computing 256 000 medium-high neurons, using one of the two logic GPUs in one K80 NVIDIA GPU. P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018)6 The forward stage eliminates the lower diagonal as follows: c0 1¼c1 b1;c0 1¼c1 bic0 i1ai fori¼2;3;...;n1 y0 1¼y1 b1;y0 1¼yiy0 i1ai bic0 i1ai fori¼2;3;...;n1 and then the backward stage recursively solves each row in reverse order: un¼y0 n;ui¼y0 ic0 iuiþ1fori¼n1;n2; :::; 1: Overall, the complexity of Thomas algorithm is optimal: 8noperations in 2n1 steps. Cyclic Reduction (CR) [7,8,20,21] is a parallel alternative to Thomas algorithm. It also consists of two phases (reduction and substitution). In each intermediate step of 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines Systems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines Systems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines Systems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines Systems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines S y stems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) 0 10 20 30 40 50 60 256 2560 25600 256000 Speedup Number of Hines S y stems Multicore Full-Inter.(one logic GPU) Full-Inter.(K80) Full-Inter.(2xK80) (a) (b) (c) (d) (e) (f) Fig. 5. Performance (speedup over sequential execution) achieved for computing multiple (256, 2560, 25 600, 256 000) neurons using different morphologies: small-low (a), small-high (b), medium-low (c), medium-high (d), big-low (e) and big-high (f). P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018) 7 the reduction phase, all even-indexed (i)equations a i x i1 +b i x i +c i x i+1 =d i are reduced. The values of a i , bi, ci and di are updated in each step according to: a0i¼ai1k1;b0 i¼bici1k1aiþ1k2 c0 i¼ciþ1k2;y0 i¼yiyi1k1yiþ1k2 k1¼ai bi1;k2¼ci biþ1: After log 2 nsteps, the system is reduced to a single equation that is solved directly. All odd-indexed unknowns xiare then solved in the substitution phase by introducing the already computed u i-1 and u i+1 values: ui¼y0 ia0 ixi1c0 ixiþ1 b0 i : Overall, the CR algorithm needs 17noperations and 2log 2 n1 steps. Figure 8a graphically illustrates its access pattern. Parallel Cyclic Reduction (PCR) [7,8,20,21]isa variant of CR, which only has substitution phase. For convenience, we consider cases where n=2 s , that involve s=log 2 nsteps. Similarly to CR, a, b, c and yare updated as follows, for j=1,2,...,sand k=2 j1 : a0i¼aiai;b0i¼biþaicikþbiaiþk c0i¼biciþ1;y0i¼biþaiyikþbiyiþk ai¼ai bi1;bi¼ci bi finally the solution is achieved as: ui¼y0 i bi: Essentially, at each reduction stage, the current system is transformed into two smaller systems and after log 2 nsteps the original system is reduced to nindependent equations. Overall, the operation count of PCR is 12nlog 2 n.Figure 8b sketches the corresponding access pattern. We should highlight that, apart from their computational complexity, these algorithms differ in their data access and synchronization patterns, which also have a strong influence on their actual performance. For instance, in the CR algorithm synchronizations are introduced at the end of each step and its corresponding memory access pattern may cause bank conflicts. PCR needs less steps and its memory access pattern is more regular [20]. In fact, hybrid combinations that try to exploit the best of each algorithm have been explored [7,8,20–23]. CR-PCR reduces the system to a certain size using the forward reduction phase of CR and then solves the reduced (intermediate) system with the PCR algorithm. Finally, it substitutes the solved unknowns back into the original system using the backward substitution phase of CR. Indeed, this is the method implemented by the gtsvStridedBatch routine into the cuSPARSE package [11], one of the implementations evaluated in this work. There are more algorithms, apart of the ones above mentioned, to deal with tridiagonal systems, such as those based on Recursive Doubling [20], among others. However, we have focused on those, which were proven to achieve a better performance and were implemented in the reference library [11]. 4.1.1 Implementation of cuThomasBatch In this section, we explore the different proposals about the CUDA thread mapping on the data layouts above 2 4 6 8 10 12 14 16 18 Mono (600-10%) Mono (600-50%) Multi (600-10%) Multi (600-50%) Speedup Fig. 6. Performance (speedup over sequential execution) achieved for computing 25 600 neurons using mono-morphologies (Mono) and multi-morphologies (Multi) of the same size and different percentages of branches (10% and 50%). ... 1 call to cuHinesBatch kernel Mono−Morphology Approach N Neurons (Hines systems) ... N Branches (Tridiagonal systems) ... N*2 Branches (Tridiagonal systems) 1 call to cuThomasBatch kernel 1 call to cuThomasBatch kernel 2 Levels Approach Multi−Morphology Fig. 7. Mono-Morphology (left) and Multi-Morphology (right) approaches. P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018)8 described on pure-tridiagonal systems. In cuThomasBatch we use a coarse-grain scheme where a set of tridiagonal systems is mapped onto a CUDA block so that each CUDA thread fully solves a system. We decided to explore this approach to avoid dealing with atomic accesses and synchronizations, as well as to be able to execute a very high number of tridiagonal systems of any size, without the limitation imposed by the parallel methods. As above presented, using the Flat data layout we cannot exploit coalescence when exploiting one thread per tridiagonal system (coarse approach); however, by interleaving (FullInterleaved data layout) the elements of the vectors, contiguous threads access to contiguous memory locations. As previously described in Section 2, this approach does not exploit efficiently the shared memory of the GPUs since the memory required by each CUDA thread becomes too large. Our GPU implementation (cuThomasBatch) is based on this approach, Thomas algorithm on Full-Interleaved data layout. On the other hand, previous studies have explored the use of the fine-grain scheme based on CR-PCR [7,8,20,21] using the Flat data layout. In this case, each tridiagonal system is distributed across the threads of a CUDA block so that the shared memory of the GPU can be used more effectively (both the matrix coefficients and the right hand side of each tridiagonal system are hold on the shared memory of the GPU). Nevertheless, computationally expensive operations, such as synchronizations and atomic accesses are necessary. Also this approach saturates the capacity of the GPU with a relatively low number of tridiagonal systems. Even when the shared memory is much faster than the global memory, it presents some important constraints to deal with. This memory is useful when the same data can be reused either by the same thread or by other thread of the same block of threads (CUDA block). Also, it is small (up to 48 KB in the architecture used) and its use hinders the exchange among blocks of threads by the CUDA scheduler to overlap accesses to global memory with computation. Our reference implementation (the gtsvStridedBatch routine into the cuSPARSE package [11]) is based on this approach, CR-PCR on Flat data layout. 4.1.2 Performance analysis To carry out the experiments, we have used one of the two logic Kepler GPUs into one K80 NVIDIA GPU. We have evaluated the performance of each of the approaches, gtsvStridedBatch and cuThomasBatch, using both, single and double precision operations. Two test cases were proposed. The first one (Figs. 9a and 10) consists of computing 256, 2560, 25 600 and 256 000 ‘‘small’’ tridiagonal systems of 64, 128, 256 and 512 elements each. Due to the memory capacity of our platform, we consider another test case (Figs. 9a and 11) for those systems with a bigger size (a higher number of elements), 1024, 2048, 4096 and 8192. In this case we could compute up to a maximum of 20 000 systems in parallel. We have considered this testbed to evaluate the scalability by increasing both, the size of the systems and the number of systems, taking into account the limitation of our platform. In particular, the size of the systems in the first test cases (64–512) can be fully executed by one CUDA block using gtsvStridedBatch. Nevertheless, those tests which need a higher size (1024 forward) must be computed following other strategies as commented before. Regarding the size of the tridiagonal systems, there is no characteristic size, as it depends on the nature of the applications, and because of that, we have considered different cases to cover all the range of possible scenarios. For the sake of numerical stability we force the tridiagonal coefficient matrix to be diagonally dominant (|b i |>|a i |+|c i |, "i=0,...,n). We initialize the matrix coefficients randomly following the previous property. Figure 9 graphically illustrates the speedup achieved by our implementation against the cuSPARSE routine. Even when interleaving the elements of the systems does not scale when computing a low number of systems (256 in Fig. 9a and 20–200 in Fig. 9b), being gtsvStridedBatch faster than our implementation, this last turns to be much faster for the rest of tests (2560–256 000 in Fig. 9a and 2000–20 000 in Fig. 9b). In most cases, independently of the size of the systems, bigger size means bigger speedup, achieving a speedup peak close to 4 in single precision and close to 3 in double precision. 2 1357 62 48 84 8642 1 345678 1 1 1 12345678 8765432 2345678 8765432 (a) (b) Fig. 8. Access pattern of the CR algorithm (a) and PCR algorithm (b). P. Valero-Lara et al.: Oil & Gas Science and Technology - Rev. IFP Energies nouvelles 73, 63 (2018) 9