# TESLA GPUs versus MPI with OpenMP for the Forward Modeling of Gravity and Gravity Gradient of Large Prisms Ensemble

Authors: Carlos Couder-Castañeda; Carlos Ortiz-Alemán; Mauricio Gabriel Orozco-del-Castillo; Mauricio Nava-Flores
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

An implementation with the CUDA technology in a single and in several graphics processing units (GPUs) is presented for the calculation of the forward modeling of gravitational fields from a tridimensional volumetric ensemble composed by unitary prisms of constant density. We compared the performance results obtained with the GPUs against a previous version coded in OpenMP with MPI, and we analyzed the results on both platforms. Today, the use of GPUs represents a breakthrough in parallel computing, which has led to the development of several applications with various applications. Nevertheless, in some applications the decomposition of the tasks is not trivial, as can be appreciated in this paper. Unlike a trivial decomposition of the domain, we proposed to decompose the problem by sets of prisms and use different memory spaces per processing CUDA core, avoiding the performance decay as a result of the constant calls to kernels functions which would be needed in a parallelization by observations points. The design and implementation created are the main contributions of this work, because the parallelization scheme implemented is not trivial. The performance results obtained are comparable to those of a small processing cluster.


---

## Body

## 1. Introduction

In recent years the number of publications about parallel computing applications using the GPUs architecture has remarkably increased. These applications represent an economic and powerful way to access high-performance computing [1, 2]. However, since the architecture of the GPU is different to that of a conventional CPU, the programming paradigm should be changed. This had led to the development of a new research field within scientific computing which explores the performance of the GPU to general purpose applications, such as acoustic simulation [3], propagation of seismic waves [4], seismic migration [5], molecular engineering [6], fluid dynamics [7], even for astrophysical simulations [8] and many other implementations. In a few words, the objective of the general purpose computing in GPU (GPGPU) is to develop new applications for those who pretend to solve problems of numerical simulation requiring as less computing time as possible.Even though the GPUs have become an accessible platform for general purpose programming, they still have some limitations and its programming entails some difficulties [2]. Compute unified device architecture (CUDA) is a set of tools that includes mainly a compiler for an extension of the C language, a set of libraries, and drivers for the specific programming of NVIDIA cards. Despite that these tools have eased the programming, it is still needed to know with precision the architecture of the card with its several memory levels to obtain the maximum performance. One of the greatest drawbacks that can occur in CUDA is the handling of critical sections or shared memory, since there are no proper instructions of exclusions, as is the case in OpenMP [13]. Therefore, when these difficulties are present in the application design, it is necessary to modify the strategy to achieve a good performance. Implementing this application in CUDA for the calculation of the direct model through a prisms ensemble represents a big challenge; this is because the memory region where the calculations are made (observation grid) is shared and therefore requires different memory allocations to avoid the data coherency problems produced when two or more processing cores access the same memory location at the same time. This drawback is not present, for example, when a model is solved through finite differences since the domain is divided between the cores and there is no overlap in the data handled by each core [9].One of the simplest parallelization options for this problem using CUDA consists of partitioning by the number of prisms and to consequently divide the domain of the observation grid between the device cores; this approach avoids the handling of shared regions. However, this design is extremely inefficient since for each prism it would be necessary to make a call to the GPU. A typical problem of fourteen million prisms would imply the execution of a kernel by the same number, and each kernel call is computationally expensive [11]. Therefore, we propose an efficient design based on the partition by groups of prisms.Additionally, we make experiments with double and single precision to calculate the errors that can introduce the single precision, and even when its use reduces the computing time from 30% to 50%, it is necessary to evaluate the effect of the introduced error by using only seven significant digits in the floating numbers and investigate if this error affects considerably the result of the modeling. This analysis is necessary since recently NVIDIA introduced TESLA K10 cards which handle single precision and are more economical than the TESLA K20 cards for double precision.


### 1.1. Related Work

Some related research works which implement an approach to calculate scalar and tensor gravity utilizing the massively parallel architecture of GPU can be found in [12], in which a parametrization based on rectilinear blocks with constant density within each block is used; however, the results show that our design yields a better performance using different memory allocations. Also a parallel program was developed to estimate the correlation imaging for gravity and gravity gradiometry data to provide a rapid approach to equivalent estimation of objective bodies with different density contrasts in the subsurface [17]; however, neither is multi-GPU implementation.


### 1.2. Paper Organization

This paper is organized as follows: in Section2 we present the characteristics of the CUDA platform and the tools we used, in Section 3 the application design is explained, in Section 4 we present some numerical experiments that were made, in Section 5 we detail the validation of the code for double as well as single precision, and finally in Section 6 we compare against a 29-nodes cluster and finally we present our conclusions.


## 1.1. Related Work

Some related research works which implement an approach to calculate scalar and tensor gravity utilizing the massively parallel architecture of GPU can be found in [12], in which a parametrization based on rectilinear blocks with constant density within each block is used; however, the results show that our design yields a better performance using different memory allocations. Also a parallel program was developed to estimate the correlation imaging for gravity and gravity gradiometry data to provide a rapid approach to equivalent estimation of objective bodies with different density contrasts in the subsurface [17]; however, neither is multi-GPU implementation.


## 1.2. Paper Organization

This paper is organized as follows: in Section2 we present the characteristics of the CUDA platform and the tools we used, in Section 3 the application design is explained, in Section 4 we present some numerical experiments that were made, in Section 5 we detail the validation of the code for double as well as single precision, and finally in Section 6 we compare against a 29-nodes cluster and finally we present our conclusions.


## 2. Architecture of the Platform

As a general-purpose architecture, CUDA includes the hardware that can have dedicated processing cards or cards which also control the visualization of the monitor and the software that includes the compiler, the card drivers, and the libraries. The programming model in CUDA consists of functions called kernels which are executed concurrently by several light threads (CUDA threads). These threads are grouped into blocks which can be of one, two, or even three dimensions. Each block can contain a maximum number of threads, defined by the architecture of the card which is being used. The blocks are executed concurrently by the stream multiprocessors (SMs) and the execution order is nondeterministic. Each SM contains a set of microprocessors which can be thought of as arithmetic logic units (ALUs) and are known as CUDA cores. The threads within a block are divided into groups of 32, called warps. A warp is executed concurrently by the CUDA cores, and the number of cores can be less than the size of the warp, as happens in the TESLA C1060 card, which has eight cores per SM but supports one warp. This configuration implies that 8 threads are executed in parallel, but 32 are processed concurrently, and this means that each core will process 4 threads previously assigned.There is an implicit synchronization between kernel calls; that is, the next kernel cannot be executed until the previous one has finalized. There are some cards of a more advanced architecture which allow the concurrent execution of kernels, but this must be specified by the programmer. The threads within the same block can be synchronized, but synchronization between blocks cannot be achieved.Understanding the different types and hierarchies of memory of the NVIDIA cards is essential to be able to take advantage of them. There are four types of memory: global, constant, texture, and shared. The global memory is analogous to the RAM memory used in a CPU. A CUDA application requires several data transfers from the global memory of the GPU to the CPU memory. The constant and the texture memories are cache memories and read-only by the SM. The content of the texture memory can be updated through special functions. The shared memory is included in each block of threads and is shared only by the threads in the block and is extremely fast in comparison with the global memory, but its deficiency lies in the fact that it is very limited and its size is defined by the architecture.The key, in general, to achieve an efficient code for GPU is to correctly handle the access times (latency) to the memory, that is, to carry out the least possible data transfers between the global memory of the GPU and the principal memory of the CPU [15], followed by few calls to the kernel functions. This implies having a great amount of blocks to process or having blocks with a great amount of threads. It is also necessary to avoid an excessive read-write access to the global memory, and preferable to use the shared memory, even though a lot of times this is not possible since a great amount of data is being handled and the shared memory becomes insufficient.


## 3. Design of the Application for GPU

The application consists of calculating the gravimetric or gradiometric response produced by a rectangular prismatic body with constant density, in reference to a set of points called observation points (see Figure1) [16]. The set of prisms is known as prisms ensemble and is not necessarily regular. An ensemble of nonregular prisms can be configured (Figure 3), with the only requisite being that they are not superimposed. Since the gravitational field complies with the superposition principle with respect to the observation, if f is the calculated response at a point (x,y), then the observed response at the point f(x,y) is given by
(1)f(x,y)=∑k=1M‍G(ρk,x,y),
where M is the total number of prisms and ρ the density of the prism.Figure 1
Construction of a prism of densities and its calculation with respect to an observation grid.It is well known that the function which can calculate the gravimetric or gradiometric contribution for a given prism and a set point can be rewritten as follows:(2)g=f(x1,y1,z1,x2,y2,z2,x,y,z,ρ),
where (x1,y1,z1) is the upper left vertex and (x2,y2,z2) the lower right vertex of the prism; (x,y,z) is the observation point and ρ the density, as shown in Figure 2.Figure 2
Illustration of the calculation of the anomaly produced by a prism with respect to an observation point.Figure 3
Illustration of an ensemble of irregular prisms.To be able to discretise the cube, we definexp, yp, and zp as the face numbers in the directions x, y, and z, respectively. If the cube is discretized in an homogeneous way, then we can define M, the number of prisms, as M=xp×yp×zp, and the consecutive numbering of the prisms would be first in x, then y and finally in z. In case the ensemble is irregular, x1, y1, z1, x2, y2, z2, and ρ should be provided through a file for each prism. We define O as the number of observation points, which is determined by O=Nx×Ny, where Nx and Ny are the number of observation points in the x and y directions, respectively. Therefore, the number of times which is required to call the function defined in (2) to calculate the anomaly produced by a component is M×O.The first step to develop a parallel program is to search the finest granularity. This is important since CUDA handles a fine granularity paradigm. In this case it can be thought of parallelizing by prisms or by observation points, (see Figure4), that is, by the number of elements in M or by the number of elements in O. One of the requisites which must be taken into consideration in the design is that it must be scalable, and hybrid systems must be considered since they are the most commonly used nowadays. However, many times design and type of architecture cannot be separated, especially when it is as specific as the CUDA architecture. Following the methodology proposed by Foster [10], it is necessary to analyze both parallelization schemes and examine which yields the best performance. Nevertheless, since M≫O, in principle the best option of partitioning is by M, as detailed in the next section, even though this design requires a greater effort for its implementation.Figure 4
Partitioning by observation points.


### 3.1. Implementation in a Single GPU

The design for a GPU will consist of generating an observation grid (memory space) for each execution thread created in the TESLA card; in other words, if we create 14 blocks, each one containing 32 execution threads, we would have generated 448 observation grids. This design obeys the fact that we will select a partition by prisms, which implies less calls to kernel functions and is therefore more efficient in terms of execution time.To give the correct dimension to our design, we analyze the simplest option of parallelization, which consists of partitioning the observation grid in the memory card for each prism.This method of parallelization is the most trivial since it is enough to simply parallelize the cycle of the observations, which can even be done by using OpenACC, simplifying even more the work and avoiding the creation of the kernel by the programmer and leaving it to the compiler. However, a big drawback of this method is the excessive number of calls to the kernel function, which decreases performance since the parallel region is created and closed on each call. Additionally, it does not represent a big challenge design-wise to the point that the scheme can be solved with a compiler which could automatically generate parallel code.On the other hand, the other parallelization option is to do it by prisms; in other words, make the threads divide the work by prisms. However, to avoid the coherence problems it is necessary to create a different space of memory for each execution thread, since it is not feasible to create just one memory space for a single observation grid, shared by all of the threads, since one of the principal problems which is not easy to handle in CUDA is the mutual exclusion of the threads in shared zones.As can be seen in Figure5, it is required to create an observation grid for each execution thread to avoid numerical consistency problems; if only one grid is occupied for all the threads, access conflicts will appear since several threads would write in the same memory location.Figure 5
Partitioning by prisms using different memory spaces.The number of created grids is equal to the number of created execution threads, and therefore the number of grids would be equal to the number of blocks times the number of threads contained in each block. This design allows the process at the same time as many prisms as threads are created, and therefore, for example, if 448 threads are created, then the same number of prisms will be processed in parallel in one kernel function execution and in the following call to the kernel function another 448 prisms will be processed and so successively until finalizing the process. In this way, the thread 1 will process the prisms set{1,449,897,1345,…}; in fact, the number of times that the kernel function is called is determined by
(3)p=⌈MT⌉,
where M is the number of prisms and T is the number of created threads, and consequently p≪M. To exemplify that the prisms partitioning is better, let us suppose that 14 blocks are created, each one containing 512 execution threads, and then a total of 7,168 observation grids are generated. If we have a problem of 200,000 prisms, the number of calls to the kernel will be ⌈200,000/7,168⌉=30, and 30 is much smaller than 200,000. Therefore, we reduced the number of calls to the kernel in a 6,666X factor with respect to the partition by observation points.Each thread will follow a scheme of processing over the prisms as follows: a threadt will process the sequence of prisms
(4)T×(n-1)+t,
where n=1,2,3,…,p.The implementation will initially consist of coding in device mode the functions which calculate the vectorial components(gx,gy,gz) and the tensorial components (gxx,gyy,gzz,gxy,gxz,gyz), following the definition of (2). The functions of type device can only be called by kernels and are executed by a single CUDA thread.To calculate any component it is necessary to allocate the memory for a tridimensional arrayG of size T×Nx×Ny. We define ID=t·x+(b·x-1)×s·x, where t·x is the thread identifier, b·x the block identifier, and s·x the block size. Notice that we write (b·x-1) because the thread and block identifiers in FORTRAN-CUDA are numbered starting from 1. We also define K as K=(T)*(I-1)+ID, where T is the number of created threads (equal to the number of observation grids), while I is the number of partition over the set of prims M in which it is working upon. The I partition of the set M is a division of M into nonoverlapping and nonempty subsets of size T that cover all of M. The subsets are collectively exhaustive and mutually exclusive with respect to the set M. The cardinality of the partition is determined by T (4) and the number of partitions is defined by p (3) and I can take the values 1,2,3,…,p. The general scheme of the computing kernel would be defined in Pseudocode 1, where Gz receives the parameters defined in (2), which are:(i)
Xa(K), Ya(K), and Za(K), the position in x, y, and z, respectively, of the upper left vertex of the prism,(ii)
Xb(K), Yb(K), and Zb(K), the position in x, y, and z, respectively, of the lower right vertex of the prism,(iii)
Xm(i), Ym(j), and Elev(i,j), the location in x, y, and z, respectively, of the observation point,(iv)
Rho(K), the density of the prism.Pseudocode 1
I
f (K<=M) then
For each
j from 1 toNy
For eachi from 1 toNx
G(Thread,i,j) =Gz(parameters) + G(Thread,i,j);
End For
End For
Else
For eachj from 1 toNy
For eachi from 1 toNx
G(Thread, i, j) = G(Thread,i,j)+0.0
End For
End For
End ifTo calculate the final result of the anomaly it is necessary to add all of the obtained results from the threads, this is, the final anomalyGf at a point (i,j) is approximated as:
(5)Gf(i,j)=∑k=1T‍G(k,i,j).To generate the reduction (the sum), we can make use of the OpenACC which automatically generates the kernel, and the structure in FORTRAN is as in Pseudocode2.Pseudocode 2
D
O k = 1,T
!$acc parallel loop present (Gd,Gshared_d)
collapse (2)
DO i = 1,Nx
DO j = 1,Ny
Gshared_D(i,j) = Gd(k,i,j) + Gshared_D(i,j);
END DO
END DO
!$acc end parallel loop
END DOIn Pseudocode2, Gshared is the bidimensional array where the reduction is made, and Gd is the tridimensional which contains the data of the anomaly generated by the different threads.


### 3.2. Multi-GPU Implementation

There are diverse possibilities to do the implementation using several GPUs [14]; one of them is to use the same libraries provided by CUDA 4.0 to the development of peer-to-peer applications to communicate the GPUs within the same workstation; another possibility would be to use OpenMP. Nevertheless, we consider that the best option in this case is to use MPI since the number of required messages is not intensive and, unlike OpenMP or the peer-to-peer connection, it allows distributing the work to several GPUs which are not even connected in the same machine. In fact, the number of messages required to be sent for a problem with q MPI processes, each one controlling a GPU, is q itself; this is because each process sends its observation grid. The estimated communication time tcomm to the complete application is
(6)tcomm=q(tstartup+tdata),
where tstartup is the necessary time to initialize the sending of the data and tdata=Nx×Ny is the size of the message, which in this case is the number of observation points.Due to the fact that the parallelization in MPI is explicit, we need to manually distribute the number of prisms through a modular expression. Let us suppose thatM is the number of prisms to calculate and that q is the number of cards that we are going to use. If every card is controlled by a process, we define the start and end of prisms to process by q as qstart and qend, respectively. We then calculate the integer s as the quotient of M and the total number of processes qn, and we determine the remainder r, both as
(7)s=Mqn,r=mod(Mqn).
Therefore
(8)qstart=q×s+1,qend=(q+1)×s.
If r≠0 and q<r, then we adjust as
(9)qstart=qstart+q,qend=qend+(q+1).
If r≠0 and q≥r, then
(10)qstart=qstart+r,qend=qend+r.
In this way we can distribute the number of prisms M over qn GPUs, in a balanced way.Once the precedent distribution is done, we can occupy the implementation of the previous section to process the subset of local prisms for each GPU. Let us suppose that we have two workstations containing 4 GPUs each one; then in each station the cards will be numbered as 0, 1, 2, 3. To correctly select a device for each processq, numbered from 0 to 7, we use the function
(11)f(q)=mod(q,nd),
where nd (constant) is the number of devices per machine (in this case 4). If q=4 then f(4)=0, and this means that the fifth process, numbered as 4, will be responsible of controlling the device 0 of the second team. It is necessary to note that this procedure works if the workload of processes is orderly distributed; in other words, in the first machine the processes 0, 1, 2, 3 are addressed and in the second 4, 5, 6, 7. If this distribution is not ordered, then the algorithm does not work properly, which is why for some MPI implementation the ordered flag can be specified as execution parameter, which orderly distributes the workload.After selecting a device per process, we proceed to distribute the prisms between the 8 processes using (7). In this way the balance for a problem of 251,946 prisms results as follows:q0={0,…,31494}
,q1={31495,…,62988},q2={62989,…,94481},q3={94482,…,125974},q4={125975,…,157467},q5={157468,…,188960},q6={188961,…,220453},q7={220454,…,251946}.After the subsets of prisms in which every card will work are defined, we apply the scheme used in Section3.1.


## 3.1. Implementation in a Single GPU

The design for a GPU will consist of generating an observation grid (memory space) for each execution thread created in the TESLA card; in other words, if we create 14 blocks, each one containing 32 execution threads, we would have generated 448 observation grids. This design obeys the fact that we will select a partition by prisms, which implies less calls to kernel functions and is therefore more efficient in terms of execution time.To give the correct dimension to our design, we analyze the simplest option of parallelization, which consists of partitioning the observation grid in the memory card for each prism.This method of parallelization is the most trivial since it is enough to simply parallelize the cycle of the observations, which can even be done by using OpenACC, simplifying even more the work and avoiding the creation of the kernel by the programmer and leaving it to the compiler. However, a big drawback of this method is the excessive number of calls to the kernel function, which decreases performance since the parallel region is created and closed on each call. Additionally, it does not represent a big challenge design-wise to the point that the scheme can be solved with a compiler which could automatically generate parallel code.On the other hand, the other parallelization option is to do it by prisms; in other words, make the threads divide the work by prisms. However, to avoid the coherence problems it is necessary to create a different space of memory for each execution thread, since it is not feasible to create just one memory space for a single observation grid, shared by all of the threads, since one of the principal problems which is not easy to handle in CUDA is the mutual exclusion of the threads in shared zones.As can be seen in Figure5, it is required to create an observation grid for each execution thread to avoid numerical consistency problems; if only one grid is occupied for all the threads, access conflicts will appear since several threads would write in the same memory location.Figure 5
Partitioning by prisms using different memory spaces.The number of created grids is equal to the number of created execution threads, and therefore the number of grids would be equal to the number of blocks times the number of threads contained in each block. This design allows the process at the same time as many prisms as threads are created, and therefore, for example, if 448 threads are created, then the same number of prisms will be processed in parallel in one kernel function execution and in the following call to the kernel function another 448 prisms will be processed and so successively until finalizing the process. In this way, the thread 1 will process the prisms set{1,449,897,1345,…}; in fact, the number of times that the kernel function is called is determined by
(3)p=⌈MT⌉,
where M is the number of prisms and T is the number of created threads, and consequently p≪M. To exemplify that the prisms partitioning is better, let us suppose that 14 blocks are created, each one containing 512 execution threads, and then a total of 7,168 observation grids are generated. If we have a problem of 200,000 prisms, the number of calls to the kernel will be ⌈200,000/7,168⌉=30, and 30 is much smaller than 200,000. Therefore, we reduced the number of calls to the kernel in a 6,666X factor with respect to the partition by observation points.Each thread will follow a scheme of processing over the prisms as follows: a threadt will process the sequence of prisms
(4)T×(n-1)+t,
where n=1,2,3,…,p.The implementation will initially consist of coding in device mode the functions which calculate the vectorial components(gx,gy,gz) and the tensorial components (gxx,gyy,gzz,gxy,gxz,gyz), following the definition of (2). The functions of type device can only be called by kernels and are executed by a single CUDA thread.To calculate any component it is necessary to allocate the memory for a tridimensional arrayG of size T×Nx×Ny. We define ID=t·x+(b·x-1)×s·x, where t·x is the thread identifier, b·x the block identifier, and s·x the block size. Notice that we write (b·x-1) because the thread and block identifiers in FORTRAN-CUDA are numbered starting from 1. We also define K as K=(T)*(I-1)+ID, where T is the number of created threads (equal to the number of observation grids), while I is the number of partition over the set of prims M in which it is working upon. The I partition of the set M is a division of M into nonoverlapping and nonempty subsets of size T that cover all of M. The subsets are collectively exhaustive and mutually exclusive with respect to the set M. The cardinality of the partition is determined by T (4) and the number of partitions is defined by p (3) and I can take the values 1,2,3,…,p. The general scheme of the computing kernel would be defined in Pseudocode 1, where Gz receives the parameters defined in (2), which are:(i)
Xa(K), Ya(K), and Za(K), the position in x, y, and z, respectively, of the upper left vertex of the prism,(ii)
Xb(K), Yb(K), and Zb(K), the position in x, y, and z, respectively, of the lower right vertex of the prism,(iii)
Xm(i), Ym(j), and Elev(i,j), the location in x, y, and z, respectively, of the observation point,(iv)
Rho(K), the density of the prism.Pseudocode 1
I
f (K<=M) then
For each
j from 1 toNy
For eachi from 1 toNx
G(Thread,i,j) =Gz(parameters) + G(Thread,i,j);
End For
End For
Else
For eachj from 1 toNy
For eachi from 1 toNx
G(Thread, i, j) = G(Thread,i,j)+0.0
End For
End For
End ifTo calculate the final result of the anomaly it is necessary to add all of the obtained results from the threads, this is, the final anomalyGf at a point (i,j) is approximated as:
(5)Gf(i,j)=∑k=1T‍G(k,i,j).To generate the reduction (the sum), we can make use of the OpenACC which automatically generates the kernel, and the structure in FORTRAN is as in Pseudocode2.Pseudocode 2
D
O k = 1,T
!$acc parallel loop present (Gd,Gshared_d)
collapse (2)
DO i = 1,Nx
DO j = 1,Ny
Gshared_D(i,j) = Gd(k,i,j) + Gshared_D(i,j);
END DO
END DO
!$acc end parallel loop
END DOIn Pseudocode2, Gshared is the bidimensional array where the reduction is made, and Gd is the tridimensional which contains the data of the anomaly generated by the different threads.


## 3.2. Multi-GPU Implementation

There are diverse possibilities to do the implementation using several GPUs [14]; one of them is to use the same libraries provided by CUDA 4.0 to the development of peer-to-peer applications to communicate the GPUs within the same workstation; another possibility would be to use OpenMP. Nevertheless, we consider that the best option in this case is to use MPI since the number of required messages is not intensive and, unlike OpenMP or the peer-to-peer connection, it allows distributing the work to several GPUs which are not even connected in the same machine. In fact, the number of messages required to be sent for a problem with q MPI processes, each one controlling a GPU, is q itself; this is because each process sends its observation grid. The estimated communication time tcomm to the complete application is
(6)tcomm=q(tstartup+tdata),
where tstartup is the necessary time to initialize the sending of the data and tdata=Nx×Ny is the size of the message, which in this case is the number of observation points.Due to the fact that the parallelization in MPI is explicit, we need to manually distribute the number of prisms through a modular expression. Let us suppose thatM is the number of prisms to calculate and that q is the number of cards that we are going to use. If every card is controlled by a process, we define the start and end of prisms to process by q as qstart and qend, respectively. We then calculate the integer s as the quotient of M and the total number of processes qn, and we determine the remainder r, both as
(7)s=Mqn,r=mod(Mqn).
Therefore
(8)qstart=q×s+1,qend=(q+1)×s.
If r≠0 and q<r, then we adjust as
(9)qstart=qstart+q,qend=qend+(q+1).
If r≠0 and q≥r, then
(10)qstart=qstart+r,qend=qend+r.
In this way we can distribute the number of prisms M over qn GPUs, in a balanced way.Once the precedent distribution is done, we can occupy the implementation of the previous section to process the subset of local prisms for each GPU. Let us suppose that we have two workstations containing 4 GPUs each one; then in each station the cards will be numbered as 0, 1, 2, 3. To correctly select a device for each processq, numbered from 0 to 7, we use the function
(11)f(q)=mod(q,nd),
where nd (constant) is the number of devices per machine (in this case 4). If q=4 then f(4)=0, and this means that the fifth process, numbered as 4, will be responsible of controlling the device 0 of the second team. It is necessary to note that this procedure works if the workload of processes is orderly distributed; in other words, in the first machine the processes 0, 1, 2, 3 are addressed and in the second 4, 5, 6, 7. If this distribution is not ordered, then the algorithm does not work properly, which is why for some MPI implementation the ordered flag can be specified as execution parameter, which orderly distributes the workload.After selecting a device per process, we proceed to distribute the prisms between the 8 processes using (7). In this way the balance for a problem of 251,946 prisms results as follows:q0={0,…,31494}
,q1={31495,…,62988},q2={62989,…,94481},q3={94482,…,125974},q4={125975,…,157467},q5={157468,…,188960},q6={188961,…,220453},q7={220454,…,251946}.After the subsets of prisms in which every card will work are defined, we apply the scheme used in Section3.1.


## 4. Performance Experiments

As experiment we use a synthetic case composed by an ensemble of700×700×50 prisms with 7 spheres of contrast of variable density (Figure 6). The spheres are conformed by 251,946 prisms and an observation grid of 150×100=15,000 points at an elevation of 100 m. Therefore, the number of calls to a function to calculate a component of the tensor or vector is 3,779,190,000, which represents a high-performance computing problem.Figure 6
Synthetic problem setup with 7 spheres of variable density contrast (not scaled). Ensemble size of 22 km× 22 km × 8 km, 251,946 conform the spheres.We performed experiments to calculate the vectorial componentsGx, Gy, and Gz, and the tensorial components Gxx, Gyy, Gzz, Gxy, Gxz, and Gyz in a single GPU, in a workstation with 4 GPUs and with two workstations containing 4 GPUs each one.The characteristics of the workstations where the experiments took place are as follows:(i)
2 Intel(R) Xeon(R) CPU X5690 @ 3.47 GHz,(ii)
6 real cores per processor,(iii)
hyperthreading technology disabled,(iv)
12 GB of RAM memory,(v)
4 TESLA Cards model C2070,(vi)
operating system Red Hat 6.3.The principal characteristics that we can highlight of the TESLA card C2070 are(i)
14 multiprocessors (SM),(ii)
32 cores per multiprocessor (448 total CUDA cores ),(iii)
6 GB of global memory DDR5,(iv)
515 theoretical GLOPS in double precision,(v)
1.03 theoretical TFLOPS in single precision,(vi)
Frequency of the CUDA cores 1.15 GHz.


### 4.1. Performance of One Single Card

The first experiment was made in a C2070 card and consists of testing with different block sizes, keeping the number of blocks fixed at 14. The size of the selected block obeys the number of SMs available in the card and varies in multiples of 32 (warp size) until 512 threads per block, which means 16 experiments.In Figure7, the obtained computing times to solve the problem of the 7 spheres in double and single precision are shown. Notice how by increasing the size of the block, the execution time decreases exponentially to a limit or an asymptotic time which in double precision is of 81 s and in single precision of 37 s. However, even when the best computing time is obtained with a block of 512 threads, this is the most memory-consuming setup. To improve the understanding of the behavior we used the speed-up as a metric, considering a block of 32 threads as the processing unit, and we label the quantity of required memory for each case.Figure 7
Comparison of the execution time using a variable block size in multiples of 32, in double and single precision.By increasing the block size in multiples of the warp, we increase what is known as multiprocessoroccupancy. In general, by increasing the occupancy we improve the use of the SM, and in this application this phenomenon is observed.The speed-up graph depicted in Figure8 shows that the computing time reduction is practically linear by increasing the occupancy for both cases, single and double precision, and later on it starts to stabilize, which means that in fact we do obtain an increase in the performance if we increase the block size in multiples of 32 threads. In the C2070 card, the maximum number of threads with their own memory space which could be created was 7,168 (14 blocks × 512 threads), with a size of 150×100=15,000 elements, which along with the number of prisms which must be previously stored in the memory of the card add up 889 MB in double precision and 471 MB in single precision. We observe that the amount of used memory in GPU increases, but the performance is improved. The number of blocks and their size—the number of created threads—will depend on the type of the card which is being used. There are several cards of medium range which do not exceed 500 MB of memory; thus so many grids cannot be created.Figure 8
Speed-up of the behavior by increasing the number of threads per block in double and single precision, with its respective quantity of global memory required.We now examine the behavior of the performance if we set the number of threads fixed at 32 per block and vary the number of blocks from 14 to 224 in multiples of 14. The results of the execution times are shown in Figure9, where it can be noted that the minimum in execution time is reached, both for single and double precision, when 112 blocks of 32 threads are created. Later on the time increases and starts decreasing again when 224 blocks are created. Comparing with the graph in Figure 7, the behavior produced by increasing the number of blocks is not as stable as increasing the number of threads per block since occupancy is not increased as only one warp per block is handled.Figure 9
Comparison of the execution time using a variable grid size in multiples of 14, in double and single precision.It is necessary to note that even though 32 blocks of 512 threads and 224 blocks of 32 threads sum up 7,168 execution threads, the second option is slower because it increases the computing time in a 42% for double precision and 54% for single precision. To simplify the writing, if we createn blocks with m threads, using double precision we write 〈n,m,d〉 and 〈n,m,s〉 in single precision.In this particular case, the creation of many blocks is not as efficient as increasing the number of threads, as shown in the speed-up graph (Figure10), considering every 14 blocks as a processing unit. In this graph a decrease of the speed-up can be clearly seen when reaching 126 blocks, but after this point the speed-up increases again until reaching at 224 blocks the same value obtained at 112 blocks.Figure 10
Speed-up of the behavior by increasing the quantity of blocks in double and single precision for a constant block size of 32 threads, with the respective quantity of required memory.


### 4.2. Performance Using Multi-GPUs

In this subsection we analyze the performance using several GPUs integrated into the same workstation and distributed in two workstations, as was mentioned in Section3.2; the choice to distribute the work between several GPUs was MPI. The best configuration found for this problem using only one C2070 card was to create 14 blocks with 512 threads, so this configuration was used. First we did experiments using four GPUs integrated into the same workstation. The results with respect to computing time are shown in Figure 11, where we considered a card as a processing unit. It can be observed that the time decrease is proportional to the number of cards and that in single precision we reduced even more the time by a factor greater than 2X. The speed-up graph shows that the decrease of time is almost linear for double and single precision (Figure 12).Figure 11
Execution time using up to four graphic cards in the same workstation in double and single precision.Figure 12
Speed-up of the behavior by increasing the number of processing cards in double and single precision. The number of blocks is 14 and the number of threads per block is 512 in every card, the optimal setup found for one single card.We then analyzed the performance using the cards in a distributed fashion. For this we used two workstations each one with four internal cards interconnected by a network of 1000 Gbytes. The results show that the communication latency is negligible since there is no perceivable overload because of the use of MPI. This is because the used MPI functions are only required at the end of the calculation to do the reductions. The configurations which are used in a distributed fashion are(1+1)    (A setup (1+1) means the use of two cards in a distributed fashion, so a configuration (m+n) implies m cards in one workstation and n in other),  (2+1), (2+2), (3+3), (4+4), and the execution times are shown in Figure 13.Figure 13
Execution time obtained using shared and distributed graphics cards, in double and single precision. We can see that MPI is not introducing overhead in the computing time.In the graph shown in Figure13 it can be observed that the execution times are practically the same for cases where the application is executed both in a local or distributed fashion. Thus, there is no difference between executing with four cards in the same workstation or with two in each workstation. With respect to the execution time, this decreases proportionally to the number of cards; however, in the cases where six and seven GPUs are used in one configuration (3+3) and (4+3), respectively, the execution times are the same since the work distribution requires the same number of executions of the kernel for each card, in this case 36. This happens because the number of prisms for the case of 7 cards is not reduced for less than 7168 per card.The behavior of the speed-up for a distributed execution is represented in Figure14. It can be observed a practically linear speed-up and in some particular cases up to a super speed-up, as is the case with the (3+3) configuration. Because of this we consider that the performance is excellent.Figure 14
Speed-up obtained using shared and distributed graphics cards, in double and single precision.


### 4.3. Comparison against a Cluster

To get a better perspective of the obtained performance with this CUDA implementation we compared it with the development made in OpenMP with MPI for cluster, and we also compared the results obtained in a cluster with the following characteristics:(i)
node: 1 Intel Xeon processors model X5550 with four physical cores per processor (2 threads per core),(ii)
29 processing nodes,(iii)
hyperthreading technology enabled,(iv)
40 GB of RAM per node,(v)
Red Hat 6.3 as operating system.The results against which this implementation is compared were obtained from a very efficient hybrid design and its implementation using MPI and OpenMP for the calculation of the gravimetry and gradiometry. The comparison of the results shows that, using 29 nodes of the cluster, it is required an execution time of 30 s; in a C2070 card using a〈14,512,d〉 configuration, the execution time is of 81 seg; four C2070 cards in the same machine communicated by MPI require 23 seg; and 8 distributed CUDA cards (4 per workstation) demand 13 s. With these times we composed the bar chart shown in Figure 15.Figure 15
Comparison of the execution times for the calculation of the tensorial components in double precision between a small cluster against a single C2070 card, four C2070 cards in the same workstation, and 8 C2070 cards distributed in two workstations.The comparison is made against 29 nodes because it is the optimal number of nodes for the distribution in a problem of 251,946 prisms with an observation grid of 15,000 points. It turns out that the cluster is 2.7X faster than a single CUDA 2070 card, but if we occupy 4 cards these are 1.3X times faster than the cluster, and we can say that there is an approximate equivalence, for this problem in particular, between 29 nodes and 4 C2070 cards, but if we occupy 8 cards distributively, these are 2.3X times faster than the cluster.


## 4.1. Performance of One Single Card

The first experiment was made in a C2070 card and consists of testing with different block sizes, keeping the number of blocks fixed at 14. The size of the selected block obeys the number of SMs available in the card and varies in multiples of 32 (warp size) until 512 threads per block, which means 16 experiments.In Figure7, the obtained computing times to solve the problem of the 7 spheres in double and single precision are shown. Notice how by increasing the size of the block, the execution time decreases exponentially to a limit or an asymptotic time which in double precision is of 81 s and in single precision of 37 s. However, even when the best computing time is obtained with a block of 512 threads, this is the most memory-consuming setup. To improve the understanding of the behavior we used the speed-up as a metric, considering a block of 32 threads as the processing unit, and we label the quantity of required memory for each case.Figure 7
Comparison of the execution time using a variable block size in multiples of 32, in double and single precision.By increasing the block size in multiples of the warp, we increase what is known as multiprocessoroccupancy. In general, by increasing the occupancy we improve the use of the SM, and in this application this phenomenon is observed.The speed-up graph depicted in Figure8 shows that the computing time reduction is practically linear by increasing the occupancy for both cases, single and double precision, and later on it starts to stabilize, which means that in fact we do obtain an increase in the performance if we increase the block size in multiples of 32 threads. In the C2070 card, the maximum number of threads with their own memory space which could be created was 7,168 (14 blocks × 512 threads), with a size of 150×100=15,000 elements, which along with the number of prisms which must be previously stored in the memory of the card add up 889 MB in double precision and 471 MB in single precision. We observe that the amount of used memory in GPU increases, but the performance is improved. The number of blocks and their size—the number of created threads—will depend on the type of the card which is being used. There are several cards of medium range which do not exceed 500 MB of memory; thus so many grids cannot be created.Figure 8
Speed-up of the behavior by increasing the number of threads per block in double and single precision, with its respective quantity of global memory required.We now examine the behavior of the performance if we set the number of threads fixed at 32 per block and vary the number of blocks from 14 to 224 in multiples of 14. The results of the execution times are shown in Figure9, where it can be noted that the minimum in execution time is reached, both for single and double precision, when 112 blocks of 32 threads are created. Later on the time increases and starts decreasing again when 224 blocks are created. Comparing with the graph in Figure 7, the behavior produced by increasing the number of blocks is not as stable as increasing the number of threads per block since occupancy is not increased as only one warp per block is handled.Figure 9
Comparison of the execution time using a variable grid size in multiples of 14, in double and single precision.It is necessary to note that even though 32 blocks of 512 threads and 224 blocks of 32 threads sum up 7,168 execution threads, the second option is slower because it increases the computing time in a 42% for double precision and 54% for single precision. To simplify the writing, if we createn blocks with m threads, using double precision we write 〈n,m,d〉 and 〈n,m,s〉 in single precision.In this particular case, the creation of many blocks is not as efficient as increasing the number of threads, as shown in the speed-up graph (Figure10), considering every 14 blocks as a processing unit. In this graph a decrease of the speed-up can be clearly seen when reaching 126 blocks, but after this point the speed-up increases again until reaching at 224 blocks the same value obtained at 112 blocks.Figure 10
Speed-up of the behavior by increasing the quantity of blocks in double and single precision for a constant block size of 32 threads, with the respective quantity of required memory.


## 4.2. Performance Using Multi-GPUs

In this subsection we analyze the performance using several GPUs integrated into the same workstation and distributed in two workstations, as was mentioned in Section3.2; the choice to distribute the work between several GPUs was MPI. The best configuration found for this problem using only one C2070 card was to create 14 blocks with 512 threads, so this configuration was used. First we did experiments using four GPUs integrated into the same workstation. The results with respect to computing time are shown in Figure 11, where we considered a card as a processing unit. It can be observed that the time decrease is proportional to the number of cards and that in single precision we reduced even more the time by a factor greater than 2X. The speed-up graph shows that the decrease of time is almost linear for double and single precision (Figure 12).Figure 11
Execution time using up to four graphic cards in the same workstation in double and single precision.Figure 12
Speed-up of the behavior by increasing the number of processing cards in double and single precision. The number of blocks is 14 and the number of threads per block is 512 in every card, the optimal setup found for one single card.We then analyzed the performance using the cards in a distributed fashion. For this we used two workstations each one with four internal cards interconnected by a network of 1000 Gbytes. The results show that the communication latency is negligible since there is no perceivable overload because of the use of MPI. This is because the used MPI functions are only required at the end of the calculation to do the reductions. The configurations which are used in a distributed fashion are(1+1)    (A setup (1+1) means the use of two cards in a distributed fashion, so a configuration (m+n) implies m cards in one workstation and n in other),  (2+1), (2+2), (3+3), (4+4), and the execution times are shown in Figure 13.Figure 13
Execution time obtained using shared and distributed graphics cards, in double and single precision. We can see that MPI is not introducing overhead in the computing time.In the graph shown in Figure13 it can be observed that the execution times are practically the same for cases where the application is executed both in a local or distributed fashion. Thus, there is no difference between executing with four cards in the same workstation or with two in each workstation. With respect to the execution time, this decreases proportionally to the number of cards; however, in the cases where six and seven GPUs are used in one configuration (3+3) and (4+3), respectively, the execution times are the same since the work distribution requires the same number of executions of the kernel for each card, in this case 36. This happens because the number of prisms for the case of 7 cards is not reduced for less than 7168 per card.The behavior of the speed-up for a distributed execution is represented in Figure14. It can be observed a practically linear speed-up and in some particular cases up to a super speed-up, as is the case with the (3+3) configuration. Because of this we consider that the performance is excellent.Figure 14
Speed-up obtained using shared and distributed graphics cards, in double and single precision.


## 4.3. Comparison against a Cluster

To get a better perspective of the obtained performance with this CUDA implementation we compared it with the development made in OpenMP with MPI for cluster, and we also compared the results obtained in a cluster with the following characteristics:(i)
node: 1 Intel Xeon processors model X5550 with four physical cores per processor (2 threads per core),(ii)
29 processing nodes,(iii)
hyperthreading technology enabled,(iv)
40 GB of RAM per node,(v)
Red Hat 6.3 as operating system.The results against which this implementation is compared were obtained from a very efficient hybrid design and its implementation using MPI and OpenMP for the calculation of the gravimetry and gradiometry. The comparison of the results shows that, using 29 nodes of the cluster, it is required an execution time of 30 s; in a C2070 card using a〈14,512,d〉 configuration, the execution time is of 81 seg; four C2070 cards in the same machine communicated by MPI require 23 seg; and 8 distributed CUDA cards (4 per workstation) demand 13 s. With these times we composed the bar chart shown in Figure 15.Figure 15
Comparison of the execution times for the calculation of the tensorial components in double precision between a small cluster against a single C2070 card, four C2070 cards in the same workstation, and 8 C2070 cards distributed in two workstations.The comparison is made against 29 nodes because it is the optimal number of nodes for the distribution in a problem of 251,946 prisms with an observation grid of 15,000 points. It turns out that the cluster is 2.7X faster than a single CUDA 2070 card, but if we occupy 4 cards these are 1.3X times faster than the cluster, and we can say that there is an approximate equivalence, for this problem in particular, between 29 nodes and 4 C2070 cards, but if we occupy 8 cards distributively, these are 2.3X times faster than the cluster.


## 5. Validation of the Numerical CUDA Code

We now proceed to verify the quality of the numerical solution produced using single card C2070 and four C2070 cards. The main objective of using the CUDA programming is to reduce the computing time; however, the quality of the numerical results must be validated and verified.To measure the error we use the formula of L2 norm, or RMS, defined as [18]
(12)e=1NxNy∑i=1Nx‍∑j=1Ny‍|gi,jp-gi,js|2,
where gi,jc is the component of the tensor calculated using the GPU, and gi,js is the component calculated serially in the CPU.In Table1 the errors of the components of the gravimetric tensor parallelly calculated with a GPU are shown, using the configuration of 〈14,512,d〉, with respect to the serial version.Table 1
Errors of the components of the gravimetric tensor, calculated with CUDA in double precision, with respect to its sequential counterpart in double precision.
Gravity components
Error L21-GPU (DP)
Error L24-GPU (DP)
G
z
2.0582
e
-
09
2.0581
e
-
09
G
x
1.7107
e
-
09
1.7106
e
-
09
G
y
1.1162
e
-
09
1.1163
e
-
09The errors in single precision for the same configuration〈14,512,s〉 are shown in Table 2.Table 2
Errors of the components of the gravimetric tensor, calculated with CUDA in single precision, with respect its sequential counterpart in double precision.
Gravity components
Error L21-GPU (SP)
Error L24-GPU (SP)
G
z
1.3080
1.3080
G
x
0.8929
0.8929
G
y
0.9275
0.9275The errors of the components of the gradient tensor of gravity calculated with CUDA in double precision, with respect to its sequential counterpart, are shown in Table3.Table 3
Errors of the components of the gradient tensor of gravity in double precision, with respect to its sequential counterpart.
Gravity gradient tensor components
Error L21-GPU (DP)
Error L24-GPU (DP)
G
z
z
4.6290
e
-
11
4.6373
e
-
11
G
x
x
7.7994
e
-
11
7.8030
e
-
11
G
y
y
8.8905
e
-
11
8.8935
e
-
11
G
x
y
5.3476
e
-
11
5.3503
e
-
11
G
x
z
9.2797
e
-
10
9.2792
e
-
10
G
y
z
2.7225
e
-
10
2.7207
e
-
10And finally in Table4 the errors of the components of the gradient tensor of gravity are shown, calculated with CUDA in single precision, with respect to its sequential counterpart.Table 4
Errors of the components of the gradient tensor of gravity in single precision, with respect to its sequential counterpart.
Gravity gradient tensor components
Error L21-GPU (SP)
Error L24-GPU (SP)
G
z
z
0.0310
0.0310
G
x
x
0.0385
0.0385
G
y
y
0.0392
0.0392
G
x
y
0.0743
0.0743
G
x
z
0.3085
0.3085
G
y
z
0.3148
0.3148It is necessary to mention that the sequential version is calculated in double precision and, as can be seen, there is practically no difference between the CUDA version in double precision and the reference solution. Nevertheless, in single precision the calculation of the vectors produces an error more significant than the calculation of the tensors. To observe how the error propagates in single precision, we show in Figure16 the calculation of Gz in double and single precision, using a single GPU.Componentz of the vector G(Gz), (a) calculated in double precision, (b) calculated in single precision. Apparently no significant differences are seen between the graphs; however there is the introduction of roughness in single precision.
(a)
(b)As can be noted, apparently the same anomaly result is reproduced both in double precision as in single, but in this last one the numerical roughness is pronounced. To examine this phenomenon in greater detail, we can zoom in Figure16 to see with more detail the introduced rugosity depicted in Figure 17.Zooming the component of the vectorGz, (a) calculated in double precision, (b) calculated in single precision. We can observe clearly the roughness introduced by the single precision.
(a)
(b)However, the numerical rugosity problem in single precision is left to the criterion of anyone interested in the result, since this precision requires practically 40% and 50% less in computing time and global memory, respectively.To note that the single precision does not always introduce pronounced numerical rugosities, we can observe the behavior of the calculation ofGzz, which does not introduce rugosity problems. Figure 18 is shows the Gzz tensor calculated in double and single precision, and Figure 19 is zoomed to focus on the detail.Componentzz of the tensor G(Gzz), (a) calculated in double precision, (b) calculated in single precision. No significant differences are noticed between the plots.
(a)
(b)Zooming the component of the tensorGzz, (a) calculated in double precision, (b) calculated in single precision. In this case the rugosity is not introduced by the single precision.
(a)
(b)Finally we mention that single precision can be used if done properly, depending on the requirements of the particular problem. The maximum absolute error found inGz with single precision is of 0.0566 mGal and for Gzz of 0.0011 Eotvos (see Figures 20 and 21). As we can see in the calculation of Gzz the introduced error is less than in Gz.Absolute error of the componentsGz calculated using four C2070 cards in (a) double precision and (b) single precision with respect to double precision CPU-cluster results.
(a)
(b)Absolute error of the componentsGzz calculated using four C2070 cards in (a) double precision and (b) single precision with respect to double precision CPU-cluster results.
(a)
(b)


## 6. Conclusions

A parallel design for the calculation of the vectorial and tensorial components of the gravity using CUDA was implemented and validated. The numerical experiments and the obtained metrics validate that the implementation is very efficient and that it also produces good results with respect to the numerical solution.We showed that selecting the simplest or most trivial parallelization technique does not necessarily leads to the best performance or the best use of the platform. In our particular case, even though the partitioning by prisms requires a greater inversion in the design and implementation, this method is the most advantageous with respect to performance because the transfers between the CPU and the GPU are minimized, moving more executable code to the GPU and also the number of calls to the kernel functions is reduced.The multi-GPU version using MPI as controller and balancer of the workload was correctly implemented since it produces practically the same results as those of the version for just one GPU. In double precision we can say that there is no difference between the calculations by the CPU and the GPU. The single precision can be used with confidence in the calculation of the tensorial components and, with appropriate considerations, in the calculation of the vectorial components as well.It was shown that for our synthetic problem, approximately 29 nodes are equivalent to 4 four C2070 cards. This obviously shows the economical benefit of using CUDA, since it is cheaper to acquire 4 graphic cards than 29 nodes of processing, and clearly the maintenance and energetic consume is considerably smaller. Nevertheless, we consider that the CUDA implementation is much more costly from the point of view of the design and the required time for its programming.We can also conclude that this implementation can serve as a design pattern to parallelize numerical schemes where the computational space cannot be disjointly divided between the processing cores, therefore minimizing the execution of the number of kernels calls.Finally we expect that GPU computing will enable us, in a near future, to optimize the numerical burden of large scale geophysical applications such as potential field modeling of impact craters [19] and multiparameter geophysical global optimization by heuristic methods [20, 21].



---

*Source: 437357-2013-11-20.xml*