# High Performance Computing of Complex Electromagnetic Algorithms Based on GPU/CPU Heterogeneous Platform and Its Applications to EM Scattering and Multilayered Medium Structure

Authors: Zhe Song; Xing Mu; Hou-Xing Zhou
Journal: International Journal of Antennas and Propagation
 (2017)

---

## Abstract

The fast and accurate numerical analysis for large-scale objects and complex structures is essential to electromagnetic simulation and design. Comparing to the exploration in EM algorithms from mathematical point of view, the computer programming realization is coordinately significant while keeping up with the development of hardware architectures. Unlike the previous parallel algorithms or those implemented by means of parallel programming on multicore CPU with OpenMP or on a cluster of computers with MPI, the new type of large-scale parallel processor based on graphics processing unit (GPU) has shown impressive ability in various scenarios of supercomputing, while its application in computational electromagnetics is especially expected. This paper introduces our recent work on high performance computing based on GPU/CPU heterogeneous platform and its application to EM scattering problems and planar multilayered medium structure, including a novel realization of OpenMP-CUDA-MLFMM, a developed ACA method and a deeply optimized CG-FFT method. With fruitful numerical examples and their obvious enhancement in efficiencies, it is convincing to keep on deeply investigating and understanding the computer hardware and their operating mechanism in the future.


---

## Body

## 1. Introduction

Demand boosting in high performance computing algorithms has been one of the most significant topics in computational electromagnetics (CEM). With the well-known merit of much fewer unknowns than finite-difference time-domain (FDTD) and finite element method (FEM), the method of moments (MoM) has been widely used in EM scattering problems of large-scale complex objects and full-wave analysis in multilayered structures during the past decades [1]. With the development in computer science and computational mathematics, many innovative algorithms have been established to accelerate the solving procedure of matrix equation. Several methods have been widely used in CEM, for example, multilevel fast multipole method (MLFMM) [2], adaptive integral method (AIM) [3, 4], fast Fourier transform based method (p-FFT, IE-FFT) [5–7], adaptive cross approximation method (ACA) [8–10], and conjugate gradient fast Fourier transform method (CG-FFT) [11]. Some of them have been successfully implemented in commercial software. Naturally, parallelization is the obvious way to make a further enhancement to efficiency [12]. In early stage, the implementation of parallel computing was mainly based on CPU platforms, such as PC-cluster with MPI and multicore CPU workstation with OpenMP. The methodology of the parallelization based on a single PC platform is the massive use of threads while sharing the same memory. The bottleneck appears immediately as the explosive growth of communication load among threads. Therefore, the theoretical speedup when using multiple processors can be predicted by the famous Amdahl’s Law [13].Thanks to the upgrading in hardware architecture of GPU, it is possible to allocate much more transistors devoted to data processing, that is, arithmetic logic unit (ALU), rather than data caching or flow control [14]. This makes the GPU most involved in highly parallel, multiprocess, and many-core computing with very high memory bandwidth. The combination of GPU and CPU can be realized by the concept of heterogeneous computing, which includes other popular hardware modules, that is, field-programmable gate array (FPGA) and digital signal processing (DSP). The programming combination between GPU and CPU can fully depend on the easy-to-use language CUDA-C/PTX. To the authors’ knowledge, the GPU/CPU heterogeneous platform is the most popular choice especially in CEM, such as the GPU-based FDTD [15], MLFMM [16–19], AIM [20], P-FFT [21], MoM [22, 23], and higher-order MoM [24]. In 2013, an impressive implementation of MLFMM by OpenMP-CUDA was realized on Fermi architecture (NVIDIA Tesla C2050), which achieved much higher performance than those before [16].In this paper, our recent works on high performance computing based on GPU/CPU heterogeneous platform are introduced with the following:(1)
A novel realization of OpenMP-CUDA-MLFMM method for EM scattering problem, in which the near-field matrix filling and the sparse matrix-vector production (MVP) are optimized, together with a warp-level parallel scheme for aggregation/disaggregation and the use of texture memory for 2D local interpolation/anterpolation(2)
A developed ACA method for EM scattering problem, in which the near-field matrix filling realizes a 100% efficiency enhancement and a thread-block-level parallel and a register-reusable schemes are applied to matrix compression and far-field MVP, respectively; also, a double-buffer technique for double precision is proposed to further enhance the efficiency of MVP computation [25](3)
A deeply optimized CG-FFT method for planar multilayered structure, in which the hardware instruction “__shufl_down” instead of shared memory was firstly used in the summation of MVP [26] and the current distribution of a large-scale flat lens was obtained efficiently.It is worth emphasizing that the domain decomposition method (DDM) is not considered in this paper, by which means all the optimization comes from arithmetic programming. Corresponding to numerical examples, the detailed parameters of GPU platform can be found in Table1, from which one can also discover the changes in computing ability brought by hardware’s upgrading. The paper is organized as follows. With a brief introduction in numerical algorithms for MLFMM, ACA, and CG-FFT, Section 2 introduces the proposed and optimized methods from the programming point of view. Section 3 deals with the analysis of fruitful numerical examples. Section 4 presents the conclusions. Some technical terms used in this paper relating to NVIDIA GPU architectures and CUDA-C language, such as thread, threadblock, grid, warp, kernel, instruction, streaming multiprocessor (SMX), and stream processor (SP), are essential for understanding and can be found in [14].Table 1
Detailed parameters of GPU platform.
Parameters
FermiGF104
KeplerGK104
KeplerGK110
Compute ability
2.1
3.0
3.5
Threads/warp
32
32
32
Max warps/multiprocessor
48
64
64
Max threads/multiprocessor
1536
2048
2048
Max thread Blocks/Multiprocessor
8
16
16
32-bit registers/multiprocessor
32768
65536
65536
Max registers/thread
63
63
255
Max threads/thread block
1024
1024
1024
Max X-grid dimension
2
16
-
1
232−1
232−1
Dynamic parallelism
×
×
√
Used in numerical algorithms
—
CG-FFT
FMM;ACA


## 2. Numerical Analysis of MLFMM, ACA, and CG-FFT

### 2.1. Principle of MLFMM and Its Optimization on GPU

As pointed out in [2], either the electric-field integral equation (EFIE) or the magnetic-field integral equation (MFIE) encounters the interior resonance at certain frequencies if exterior medium is lossless. The combined-field integral equation (CFIE) is always used to eliminate the problem, in which the matrix element can be calculated by(1)Zij=αηf⃑i,Kf⃑j+1-αf⃑i,Lf⃑j,where the EFIE and MFIE are presented as(2)f⃑i,Lf⃑j=jkη∫Sif⃑ir⃑i·ds∫SjI̿+1k2∇∇Gr⃑i,r⃑j′·f⃑jr⃑j′ds′,f⃑i,Kf⃑j=-∫Sif⃑ir⃑i×n^i·ds∫Sj∇Gr⃑i,r⃑j′×f⃑jr⃑j′ds′.By using the addition theorem of Green’s function, formula (1) can be rewritten in the form named as “aggregation–translation–disaggregation” for EFIE and MFIE:(3)f⃑i,Df⃑j=∫Ed2k^V⃑ml′,jS,Pk^·αmlml′k^·V⃑ml′,iF,Pk^,V⃑ml′,jS,Dk^=∫Sjds′I̿-k^k^·f⃑jr⃑j′e-jk⃑·r⃑ml′j,αmlml′k^=jkη4πTBk^,r⃑mlml′,V⃑ml,iF,Dk^=∫SidsI̿-k^k^·f⃑ir⃑ie-jk⃑·r⃑iml,TBk^,r⃑≐k4π∑l=0B-jl+1·2l+1hl2krPlk^·r⃑,where D=K or L, hl2·, and Pl· stand for thel-order Hankel function of the second kind and thel-order Legendre function, respectively. TB is the truncation length and related to the diameter of the cubic in lth layer and the numerical accuracy.Since the generation of near-field matrix is a highly intensive work in computation of MLFMM, the implementation of GPU or multi-GPU can significantly speedup the procedure even with RWG-to-RWG rather than triangle-to-triangle scheme. All the RWG functions are stored as two copies, namely, testing chain and basis chain, and each thread on GPU deals with one or more nonzero matrix elements by accessing corresponding data from the chains. This scheme can easily extend to multi-GPU case by splitting the whole task into subtasks. For example, if there are 2 GPUs, namely, GPU-0 and GPU-1, the chain can be split asfiti=0N/2-1and fiti=N/2N, which makes each GPU has a N/2×N sparse matrix to deal with. Considering the asynchronous launching scheme of kernel functions in CUDA, the control returns to CPU after launch immediately, which means a duty waste for CPU when a thread is working on GPU. Therefore, it is necessary to establish a heterogeneous scheme for coordination between CPU and GPU. For this reason, each near matrix is further divided into two parts withNG andNC rows, respectively, as shown in Figure 1. All the tasks can be executed concurrently while satisfying the asynchronous mechanism of GPU/CPU heterogeneous platform. This kind of subdivisions can accelerate the sparse matrix-vector productions (SpMVP). The determination ofNG andNC cannot be easily derived from mathematics; however, the empirical value can be obtained from numerical simulation and curve fitting, in which the optimum solution can be refined by parameter scanning, as shown in Figure 2. Considering the coexistence of multicore CPU and multi-GPU,nC andnG stand for the number of cores in CPU and number of GPU, respectively. The sampling points were obtained from sphere scatterings with radius changing from one to six wavelengths. As a result, for the case ofnC = 4 andnG = 2, the optimum solution ofNC/(N/2) equals 0.016 and 0.083 for single and double precision, respectively.Figure 1
The subdivisions of the near matrix and the task assignments.Figure 2
N
C as a function ofN/nG with nG=2.
(a)
Single precision
(b)
Double precisionThe implementation of aggregation and disaggregation at finest level on GPU was proposed by means of allocating a thread to each spectrum point [16]. To increase the utilization of GPU on Kepler architecture (GK110) further, which has a maximum value of 32 threads in one warp, a smart scheme is designed with two steps. Firstly, all the cubes are grouped into each GPU evenly. Secondly, one warp is assigned to a certain cube, in which one thread is assigned to one or more RWG functions immediately. Taking the advantage of the hardware architecture in Kepler GK110, the shuffle instruction “__shulf()” can efficiently calculate the reductive summation in register rather than shared memory [14]. Meanwhile, all the 32 threads in a warp can read the data from constant memory according to a certain spectrum point through the read-only cache by using the instruction “__ldg()” or “__restrict__”.Unlike the strategy of thread-based task assignment proposed for aggregation and disaggregation at coarser level [16], we make a further step in data storage by using the four times larger texture memory in Kepler than Fermi. Since the local inter-/anterpolation accesses neighboring data frequently, it had better store data in a certain pattern form looking like the geometric topology in texture memory. As is pointed out in [27], the texture cache is optimized for 2D spatial locality. Therefore, the threads in one warp reading close-set texture addresses can acquire hardware acceleration and achieve the best performance. With texture memory, the process from coarser level towards finest level can also be executed efficiently. In detail, after all the cubes at child level are grouped evenly in every device (i.e., GPU), a thread will be assigned to give a cube for a spectrum point at this level. By fetching the data from global memory, the information can be organized into a 3D CUDA array that is bound to the layer of texture memory. Then, the cube receives data disaggregated from its parent level by applying local anterpolation within texture memory. All this kind of tasks for cubes at child level can be executed concurrently and all the layers of texture memory assigned can be flushed for the subsequent requests.


### 2.2. Principle of ACA and Its Optimization on GPU

The ACA is a purely algebraic algorithm [8, 9], aiming at independence from the integral kernel, that is, Green’s function, which makes it different from MLFMM. The main idea of ACA is based on the compressing of a rank deficient system, which corresponds to the submatrix representing interactions between two well-separated groups of RWGs. Therefore, the MVP of this kind of submatrices can be efficiently calculated by using the factorial forms [9]. In ACA, the MVP is divided into near- and far-field interactions, which leads to the same data structure as that in MLFMM, that is, the Octree structure, and forms a hierarchical representation of far-field submatrices. Hereby we use Zm,nl to denote the m×n submatrix for interaction between groupm and groupn at the lth-level with indicating a well-separated situation. Then, it can be approximated as(4)Zm,nl≈Z~m,nl=Um,klVk,nl,(5)k≈rankUm,nl<minm,nwhich satisfies(6)Zm,nl-Um,klVk,nlF<εACAZm,nlF,where ·F is the Frobenius norm of a matrix and εACA is a given threshold. The workflow of ACA can be depicted as in the following steps with initial value ofk = 1 andik = 1.Step 1.
Generatevk as the ith row of Zm,nl and update it.(7)vkj⟸vkj-∑i=1k-1uikvij.Step 2.
Ifvkj=0,∀j, then ik⇐ik+1 and go back to Step 1; otherwise, normalize the vector with its maximum element,(8)vksk=maxj⁡vkj,vkj⟸vkjvksk,∀j.Step 3.
Generateuk as theskth row of Zm,nl and update it,(9)ukj⟸ukj-∑i=1k-1uijvisk.Step 4.
GenerateUm,kl and Vk,nl, where(10)Um,kl=u1,…,uk,Vk,nl=v1⋮vk.Step 5.
The feedback system will keep updating byk⇐k+1 and ik⇐ik+1 until the condition in (5) is satisfied.Obviously, it is convenient to employ both threadgrid and threadblock to take charge of one dimension. The total number of threads in each block is selected as 256, by means of 8 warps. The RWG functions, no matter what it is seen as, basis or test function, can be stored in shared memory during the generation of factorial form of a submatrix. Meanwhile, the submatricesU andV are stored in column- and row-major order, respectively. For the task assignment, 2l warps perform the ACA calculation in one submatrix at the (L-l)th level. However, this arrangement has no worry about insufficient blocks or logic error forl = 0, 1, 2, 3 untill = 4, because there is no synchronization mechanism between different threadblocks. In this paper, we realize one threadblock to one submatrix even at (l ≥ 4) by launching a child CUDA kernel function during the generation of new row or column vector [25]. The number of threads enabled in this so-called child kernel is equal to the length of the row or column vector and they calculate only one entry of the vector. It is worth emphasizing that an explicit synchronization between parent and child kernels should be adjusted in advance, which ensures that all the tasks assigned to child kernel are ultimately completed before logic and data operates on parent kernel. This treatment can be realized by the specific “dynamic parallelism” [28] in Kepler GK110 (Table 1) and a vivid flowchart is shown in Figure 3.Figure 3
Flowchart of dynamic parallelism.After the ACA operation, the MoM matrix is converted into anH-matrix that greatly reduces the storage and speeds up the MVPs. The upcoming procedure is to calculate massive MVPs, which is still time consuming even though theH-matrix is used. In this paper, the batched matrix-vector productions (bMVP) are considered to deal with the massive computation [25]. The main idea of bMVP is to apply the 2D form of threadblock with dim3(pth, gridDimy) to calculate all the MVPs. The column inU and the row inV are performed concurrently, such that all segments (column or row) are with the same size except the last part, as shown in Figure 4. The result of y=Vk,nlix is stored in the shared memory as an intermediate buffer and then the result of z=Um,kliy is accumulated into the result by reductive summation in shared memory. If we lets be equal to 256, there will be 64 threads in a group with the same ID iny-direction to do multiplication and the vector x can be reused in registers. Therefore, each threadblock can deal with 8 rows ofV matrix and at most 32 threadblocks (gridDimy = 32) are needed to calculate the MVP.Figure 4
The segmentation scheme.


### 2.3. Principle of CG-FFT, Multilayered Green’s Functions, and Its Optimization on GPU

The CG-FFT [11] is one of the most classical methods with iterative procedure and many applications have been found in different areas such as algebra, signal processing, computational geometry, and CEM, especially in MoM based algorithms [29]. The pseudocode of traditional CG-FFT can be depicted asx=x0, initial valuer0=b-Ax0,p0=s0=AHr0,γ0=s022fork=0,1,2, until convergenceqk=Apk, MVPαk=γk/qk22,xk+1=xk+αkpk,rk+1=rk-αkqk,sk+1=AHrk+1, MVPγk+1=sk+122,βk=γk+1/γk,pk+1=sk+1+βkpk,where AH is the conjugate transpose of matrixA andk stands for the iteration round. Obviously, the MVP operation dominates the most computation source and needs to be optimized. As we mentioned before, after a multiprocessor is assigned one or more threadblocks to execute, it will divide them into warps and each warp will be scheduled by a warp scheduler. The full efficiency can be realized when all 32 threads of a warp agree on their execution path. Considering the Kepler GK110 platform, the maximum threadblocks in one stream multiprocessor are 16 and each SMX contains 192 ALUs. An optimized parallelization can be realized by the so-called “single instruction multiple thread (SIMT)” [14]. The flowchart of this optimized CG-FFT on GPU is depicted in Figure 5, in which the “RHT” stands for right-hand term in MoM matrix equation.Figure 5
Flowchart of matrix partition for threads assignment.As the matrix is partitioned intoN submatrices with at least (N-1) of them having the same size, each submatrix can also be partitioned into pieces and each piece occupies a certain thread. Just like the realization in OpenMP-CUDA-MLFMM, the hardware instruction “__shufl()” is used instead of the shared memory to do reductive summation. This makes the performance even more efficient than the famous CUBLAS library. That is because the shuffle instruction is more efficient and reduces the usage of shared memory with increasing the occupancy.A brief review of dyadic Green’s functions in planar multilayered medium structure is also necessary in this part. The most important one is the mixed potential integral equation (MPIE), based on which the MoM can be established [30–35].(11)E⃑=-jωμ0G̿A;J⃑+1jωε0∇GΦ;∇′·J⃑,where the dyadic Green’s functions can be presented as (12)G̿A=GxxA0GxzA0GxxAGyzAGzxAGzyAGzzA.One of the most efficient ways to evaluate the Green’s functions through Sommerfeld Integral (SI) from the spectral to spatial domain is the discrete complex image method (DCIM). According to the analysis in [30], the DCIM is suitable for near-field region (ρ ≤ 0.05λ) by calculating three terms, namely, the quasi-static term (GQ-S), the surface wave term (GSWP), and the complex image term (GCIM), as shown below.(13a)Gρ=A4πGQ-S+GSWP+GCIM,(13b)GQ-S=Gρ→∞·e-jk0ρρ,(13c)GSWP=-2πj∑iNSWPResiH02kρiρkρi,(13d)GCIM=∑i=1Ncaie-jk0riri,where H02 stand for the zero-order Hankel function of the second kind and Resi stands for the residues of surface wave poles on the complex spectrum plane. For far-field region (ρ ≥ 0.05λ), both the surface and leaky wave should be considered, together with a branch cut integral contribution, as shown below.(14a)∫SIP∙dkρ=∫Γ∙dkρ-2πj·∑iNSWPRkρi+∑jNLWPRkρj,(14b)∫Γ∙dkρ=∫ΓG~+kρ+G~-kρ·Hn2kρρkρdkρ,where SIP stands for the Sommerfeld integral path, NSWP and NLWP stand for the number of extracted surface and leaky wave poles, respectively, Г stands for the integral path along the branch cut, and G~+ and G~- are the spectral Green’s functions along the top and bottom Riemann sheets of the branch cut, respectively [30]. With accurate evaluation of all the components in dyadic Green’s functions, the MPIE can be solved by MoM with RWG functions, in which the matrix element is calculated as(15)Zmn=k02∑n=1NIn∫Tm∫Tnf⃑nr⃑′·G̿Ar⃑∣r⃑′·f⃑mr⃑ds′ds+∑n=1NIn·∫Tm∇∫TnGΦr⃑∣r⃑′∇S′·f⃑nr⃑′ds′·f⃑mr⃑ds,where f⃑n and f⃑m stand for the basis and testing RWG functions, respectively. The right-hand term of the MoM matrix equation will be different from radiating, transmitting to scattering cases.


## 2.1. Principle of MLFMM and Its Optimization on GPU

As pointed out in [2], either the electric-field integral equation (EFIE) or the magnetic-field integral equation (MFIE) encounters the interior resonance at certain frequencies if exterior medium is lossless. The combined-field integral equation (CFIE) is always used to eliminate the problem, in which the matrix element can be calculated by(1)Zij=αηf⃑i,Kf⃑j+1-αf⃑i,Lf⃑j,where the EFIE and MFIE are presented as(2)f⃑i,Lf⃑j=jkη∫Sif⃑ir⃑i·ds∫SjI̿+1k2∇∇Gr⃑i,r⃑j′·f⃑jr⃑j′ds′,f⃑i,Kf⃑j=-∫Sif⃑ir⃑i×n^i·ds∫Sj∇Gr⃑i,r⃑j′×f⃑jr⃑j′ds′.By using the addition theorem of Green’s function, formula (1) can be rewritten in the form named as “aggregation–translation–disaggregation” for EFIE and MFIE:(3)f⃑i,Df⃑j=∫Ed2k^V⃑ml′,jS,Pk^·αmlml′k^·V⃑ml′,iF,Pk^,V⃑ml′,jS,Dk^=∫Sjds′I̿-k^k^·f⃑jr⃑j′e-jk⃑·r⃑ml′j,αmlml′k^=jkη4πTBk^,r⃑mlml′,V⃑ml,iF,Dk^=∫SidsI̿-k^k^·f⃑ir⃑ie-jk⃑·r⃑iml,TBk^,r⃑≐k4π∑l=0B-jl+1·2l+1hl2krPlk^·r⃑,where D=K or L, hl2·, and Pl· stand for thel-order Hankel function of the second kind and thel-order Legendre function, respectively. TB is the truncation length and related to the diameter of the cubic in lth layer and the numerical accuracy.Since the generation of near-field matrix is a highly intensive work in computation of MLFMM, the implementation of GPU or multi-GPU can significantly speedup the procedure even with RWG-to-RWG rather than triangle-to-triangle scheme. All the RWG functions are stored as two copies, namely, testing chain and basis chain, and each thread on GPU deals with one or more nonzero matrix elements by accessing corresponding data from the chains. This scheme can easily extend to multi-GPU case by splitting the whole task into subtasks. For example, if there are 2 GPUs, namely, GPU-0 and GPU-1, the chain can be split asfiti=0N/2-1and fiti=N/2N, which makes each GPU has a N/2×N sparse matrix to deal with. Considering the asynchronous launching scheme of kernel functions in CUDA, the control returns to CPU after launch immediately, which means a duty waste for CPU when a thread is working on GPU. Therefore, it is necessary to establish a heterogeneous scheme for coordination between CPU and GPU. For this reason, each near matrix is further divided into two parts withNG andNC rows, respectively, as shown in Figure 1. All the tasks can be executed concurrently while satisfying the asynchronous mechanism of GPU/CPU heterogeneous platform. This kind of subdivisions can accelerate the sparse matrix-vector productions (SpMVP). The determination ofNG andNC cannot be easily derived from mathematics; however, the empirical value can be obtained from numerical simulation and curve fitting, in which the optimum solution can be refined by parameter scanning, as shown in Figure 2. Considering the coexistence of multicore CPU and multi-GPU,nC andnG stand for the number of cores in CPU and number of GPU, respectively. The sampling points were obtained from sphere scatterings with radius changing from one to six wavelengths. As a result, for the case ofnC = 4 andnG = 2, the optimum solution ofNC/(N/2) equals 0.016 and 0.083 for single and double precision, respectively.Figure 1
The subdivisions of the near matrix and the task assignments.Figure 2
N
C as a function ofN/nG with nG=2.
(a)
Single precision
(b)
Double precisionThe implementation of aggregation and disaggregation at finest level on GPU was proposed by means of allocating a thread to each spectrum point [16]. To increase the utilization of GPU on Kepler architecture (GK110) further, which has a maximum value of 32 threads in one warp, a smart scheme is designed with two steps. Firstly, all the cubes are grouped into each GPU evenly. Secondly, one warp is assigned to a certain cube, in which one thread is assigned to one or more RWG functions immediately. Taking the advantage of the hardware architecture in Kepler GK110, the shuffle instruction “__shulf()” can efficiently calculate the reductive summation in register rather than shared memory [14]. Meanwhile, all the 32 threads in a warp can read the data from constant memory according to a certain spectrum point through the read-only cache by using the instruction “__ldg()” or “__restrict__”.Unlike the strategy of thread-based task assignment proposed for aggregation and disaggregation at coarser level [16], we make a further step in data storage by using the four times larger texture memory in Kepler than Fermi. Since the local inter-/anterpolation accesses neighboring data frequently, it had better store data in a certain pattern form looking like the geometric topology in texture memory. As is pointed out in [27], the texture cache is optimized for 2D spatial locality. Therefore, the threads in one warp reading close-set texture addresses can acquire hardware acceleration and achieve the best performance. With texture memory, the process from coarser level towards finest level can also be executed efficiently. In detail, after all the cubes at child level are grouped evenly in every device (i.e., GPU), a thread will be assigned to give a cube for a spectrum point at this level. By fetching the data from global memory, the information can be organized into a 3D CUDA array that is bound to the layer of texture memory. Then, the cube receives data disaggregated from its parent level by applying local anterpolation within texture memory. All this kind of tasks for cubes at child level can be executed concurrently and all the layers of texture memory assigned can be flushed for the subsequent requests.


## 2.2. Principle of ACA and Its Optimization on GPU

The ACA is a purely algebraic algorithm [8, 9], aiming at independence from the integral kernel, that is, Green’s function, which makes it different from MLFMM. The main idea of ACA is based on the compressing of a rank deficient system, which corresponds to the submatrix representing interactions between two well-separated groups of RWGs. Therefore, the MVP of this kind of submatrices can be efficiently calculated by using the factorial forms [9]. In ACA, the MVP is divided into near- and far-field interactions, which leads to the same data structure as that in MLFMM, that is, the Octree structure, and forms a hierarchical representation of far-field submatrices. Hereby we use Zm,nl to denote the m×n submatrix for interaction between groupm and groupn at the lth-level with indicating a well-separated situation. Then, it can be approximated as(4)Zm,nl≈Z~m,nl=Um,klVk,nl,(5)k≈rankUm,nl<minm,nwhich satisfies(6)Zm,nl-Um,klVk,nlF<εACAZm,nlF,where ·F is the Frobenius norm of a matrix and εACA is a given threshold. The workflow of ACA can be depicted as in the following steps with initial value ofk = 1 andik = 1.Step 1.
Generatevk as the ith row of Zm,nl and update it.(7)vkj⟸vkj-∑i=1k-1uikvij.Step 2.
Ifvkj=0,∀j, then ik⇐ik+1 and go back to Step 1; otherwise, normalize the vector with its maximum element,(8)vksk=maxj⁡vkj,vkj⟸vkjvksk,∀j.Step 3.
Generateuk as theskth row of Zm,nl and update it,(9)ukj⟸ukj-∑i=1k-1uijvisk.Step 4.
GenerateUm,kl and Vk,nl, where(10)Um,kl=u1,…,uk,Vk,nl=v1⋮vk.Step 5.
The feedback system will keep updating byk⇐k+1 and ik⇐ik+1 until the condition in (5) is satisfied.Obviously, it is convenient to employ both threadgrid and threadblock to take charge of one dimension. The total number of threads in each block is selected as 256, by means of 8 warps. The RWG functions, no matter what it is seen as, basis or test function, can be stored in shared memory during the generation of factorial form of a submatrix. Meanwhile, the submatricesU andV are stored in column- and row-major order, respectively. For the task assignment, 2l warps perform the ACA calculation in one submatrix at the (L-l)th level. However, this arrangement has no worry about insufficient blocks or logic error forl = 0, 1, 2, 3 untill = 4, because there is no synchronization mechanism between different threadblocks. In this paper, we realize one threadblock to one submatrix even at (l ≥ 4) by launching a child CUDA kernel function during the generation of new row or column vector [25]. The number of threads enabled in this so-called child kernel is equal to the length of the row or column vector and they calculate only one entry of the vector. It is worth emphasizing that an explicit synchronization between parent and child kernels should be adjusted in advance, which ensures that all the tasks assigned to child kernel are ultimately completed before logic and data operates on parent kernel. This treatment can be realized by the specific “dynamic parallelism” [28] in Kepler GK110 (Table 1) and a vivid flowchart is shown in Figure 3.Figure 3
Flowchart of dynamic parallelism.After the ACA operation, the MoM matrix is converted into anH-matrix that greatly reduces the storage and speeds up the MVPs. The upcoming procedure is to calculate massive MVPs, which is still time consuming even though theH-matrix is used. In this paper, the batched matrix-vector productions (bMVP) are considered to deal with the massive computation [25]. The main idea of bMVP is to apply the 2D form of threadblock with dim3(pth, gridDimy) to calculate all the MVPs. The column inU and the row inV are performed concurrently, such that all segments (column or row) are with the same size except the last part, as shown in Figure 4. The result of y=Vk,nlix is stored in the shared memory as an intermediate buffer and then the result of z=Um,kliy is accumulated into the result by reductive summation in shared memory. If we lets be equal to 256, there will be 64 threads in a group with the same ID iny-direction to do multiplication and the vector x can be reused in registers. Therefore, each threadblock can deal with 8 rows ofV matrix and at most 32 threadblocks (gridDimy = 32) are needed to calculate the MVP.Figure 4
The segmentation scheme.


## 2.3. Principle of CG-FFT, Multilayered Green’s Functions, and Its Optimization on GPU

The CG-FFT [11] is one of the most classical methods with iterative procedure and many applications have been found in different areas such as algebra, signal processing, computational geometry, and CEM, especially in MoM based algorithms [29]. The pseudocode of traditional CG-FFT can be depicted asx=x0, initial valuer0=b-Ax0,p0=s0=AHr0,γ0=s022fork=0,1,2, until convergenceqk=Apk, MVPαk=γk/qk22,xk+1=xk+αkpk,rk+1=rk-αkqk,sk+1=AHrk+1, MVPγk+1=sk+122,βk=γk+1/γk,pk+1=sk+1+βkpk,where AH is the conjugate transpose of matrixA andk stands for the iteration round. Obviously, the MVP operation dominates the most computation source and needs to be optimized. As we mentioned before, after a multiprocessor is assigned one or more threadblocks to execute, it will divide them into warps and each warp will be scheduled by a warp scheduler. The full efficiency can be realized when all 32 threads of a warp agree on their execution path. Considering the Kepler GK110 platform, the maximum threadblocks in one stream multiprocessor are 16 and each SMX contains 192 ALUs. An optimized parallelization can be realized by the so-called “single instruction multiple thread (SIMT)” [14]. The flowchart of this optimized CG-FFT on GPU is depicted in Figure 5, in which the “RHT” stands for right-hand term in MoM matrix equation.Figure 5
Flowchart of matrix partition for threads assignment.As the matrix is partitioned intoN submatrices with at least (N-1) of them having the same size, each submatrix can also be partitioned into pieces and each piece occupies a certain thread. Just like the realization in OpenMP-CUDA-MLFMM, the hardware instruction “__shufl()” is used instead of the shared memory to do reductive summation. This makes the performance even more efficient than the famous CUBLAS library. That is because the shuffle instruction is more efficient and reduces the usage of shared memory with increasing the occupancy.A brief review of dyadic Green’s functions in planar multilayered medium structure is also necessary in this part. The most important one is the mixed potential integral equation (MPIE), based on which the MoM can be established [30–35].(11)E⃑=-jωμ0G̿A;J⃑+1jωε0∇GΦ;∇′·J⃑,where the dyadic Green’s functions can be presented as (12)G̿A=GxxA0GxzA0GxxAGyzAGzxAGzyAGzzA.One of the most efficient ways to evaluate the Green’s functions through Sommerfeld Integral (SI) from the spectral to spatial domain is the discrete complex image method (DCIM). According to the analysis in [30], the DCIM is suitable for near-field region (ρ ≤ 0.05λ) by calculating three terms, namely, the quasi-static term (GQ-S), the surface wave term (GSWP), and the complex image term (GCIM), as shown below.(13a)Gρ=A4πGQ-S+GSWP+GCIM,(13b)GQ-S=Gρ→∞·e-jk0ρρ,(13c)GSWP=-2πj∑iNSWPResiH02kρiρkρi,(13d)GCIM=∑i=1Ncaie-jk0riri,where H02 stand for the zero-order Hankel function of the second kind and Resi stands for the residues of surface wave poles on the complex spectrum plane. For far-field region (ρ ≥ 0.05λ), both the surface and leaky wave should be considered, together with a branch cut integral contribution, as shown below.(14a)∫SIP∙dkρ=∫Γ∙dkρ-2πj·∑iNSWPRkρi+∑jNLWPRkρj,(14b)∫Γ∙dkρ=∫ΓG~+kρ+G~-kρ·Hn2kρρkρdkρ,where SIP stands for the Sommerfeld integral path, NSWP and NLWP stand for the number of extracted surface and leaky wave poles, respectively, Г stands for the integral path along the branch cut, and G~+ and G~- are the spectral Green’s functions along the top and bottom Riemann sheets of the branch cut, respectively [30]. With accurate evaluation of all the components in dyadic Green’s functions, the MPIE can be solved by MoM with RWG functions, in which the matrix element is calculated as(15)Zmn=k02∑n=1NIn∫Tm∫Tnf⃑nr⃑′·G̿Ar⃑∣r⃑′·f⃑mr⃑ds′ds+∑n=1NIn·∫Tm∇∫TnGΦr⃑∣r⃑′∇S′·f⃑nr⃑′ds′·f⃑mr⃑ds,where f⃑n and f⃑m stand for the basis and testing RWG functions, respectively. The right-hand term of the MoM matrix equation will be different from radiating, transmitting to scattering cases.


## 3. Numerical Examples

To verify the efficiency enhancement of the proposed optimization based on GPU/CPU heterogeneous platform, this section will demonstrate fruitful numerical examples, respectively, corresponding to the OpenMP-CUDA-MLFMM, the GPU-based ACA, and the deeply optimized CG-FFT. All the programming codes are compiled in CUDA-C environment, and the hardware configuration can be found in Table2. Similar to Section 2, the numerical results are also separated into 3 groups. The correctness of proposed methods with MLFMM and ACA have been verified by the very good agreement between bistatic RCS and standard Mie solutions of PEC spheres, which will not be demonstrated in this section.Table 2
Configurations of GPU/CPU platform of numerical examples.
Parameters
MLFMM
ACA
CG-FFT
CPU
Intel
Intel
Intel Xeon
Processor #
i5 4570
i5 4570
E5420
# of cores
4
4
4
Base frequency (GHz)
3.20
3.20
2.50
GPU
NVIDIA
NVIDIA
NVIDIA
Architecture
GK110
GK110
GK104
CUDA cores
2304
2304
960
Base clock (MHz)
863
863
980
Memory Config (MB)
3072
3072
2048
Memory BW (GB/s)
288
288
144
CUDA version
6.0
7.0 RC
4.2
# of GPU
2
2
1


### 3.1. Numerical Results of OpenMP-CUDA-MLFMM

As a benchmark result to verify the robustness of algorithms, the scattering from the PEC NASA Almond at 7 GHz for VV polarization and that from the PEC NASA Ogive at 9 GHz for HH polarization are considered. These two PEC surfaces are discretized into refined triangles with edge length of about 0.1λ, which produce 8416 and 15516 RWG unknowns, respectively. The monostatic RCS curves with the incident angle as argument are obtained by the proposed method, as shown in Figure 6. From the comparison with measurement results, very good agreement can be found. Note that the interpretations labeled as “mid VV,” “high VV,” and “VV FERM” refer to Figure 5 in [36], and those labeled as “HH,” “VV,” “HH Cicero,” and “VV Cicero” refer to Figure 9 in [36]. To demonstrate the efficiency performance of the proposed OpenMP-CUDA-MLFMM, we need to define some variables to record time consumption in different procedures, as shown in Table 3.Table 3
Variables definition of time consumption.
Variables
Description
t
M
Time for near matrix generation
t
S
Time for corresponding SpMVP
t
A
F
Time for aggregation at finest level
t
D
F
Time for disaggregation at finest level
t
A
S
Time for aggregation towards coarsest level
t
D
S
Time for disaggregation towards finest levelFigure 6
(a) The monostatic RCS at 7 GHz for VV polarization. (b) The monostatic RCS at 9 GHz for HH polarization.
(a)
(b)By sampling on different number of unknowns, an approximate curve of(tM;tS) and (tAF+tDF;tAS+tDS) can be fitted in a broad range, as shown in Figure 7. In can be clearly seen that the proposed method has higher efficiency than that reported in [16] based on the same Kepler GK110 architecture. It is worth pointing out that the translation is not included because both the proposed method and that in [16] have almost the same performance in this procedure.Figure 7
t
M and tS as function ofN with (a) single and (b) double precision; tAF+tDF and tAS+tDS as function ofN with (c) single and (d) double precision.
(a)
(b)
(c)
(d)


### 3.2. Numerical Results of GPU-Based ACA

To verify the robustness of the proposed GPU-based ACA, the bistatic RCS from a rotor with four blades and a depressed cylindrical object are considered, as shown in Figure8. Very good agreement can be seen between the proposed method and the traditional MoM on CPU. These two PEC surfaces are discretized into refined triangles with edge length of about 0.1λ, which produce 6636 and 4817 RWG unknowns, respectively.Figure 8
Bistatic RCS of (a) 4 blades rotor and (b) depressed cylindrical object.
(a)
(b)To demonstrate the efficiency performance of the proposed GPU-based ACA, the scattering from a PEC sphere with radius ofnλ (n= 1, 2, 3, and 4) is considered. The number of generated RWG unknowns is 4774, 18468, 40950, and 72120, respectively. In this paper, the threshold in formula (6) is selected as εACA=10-3. The efficiency comparison between GPU and CPU platform in both the bMVP operations and the whole ACA scheme with single and double precision are shown in Figure 9. It can be seen that the speedup ratio of matrix compression by using ACA can achieve about 50~100 and 25~50 for single and double precision, respectively. Meanwhile, the speedup ratio of bMVP can achieve about  10~30 and 6~17 for single and double precision, respectively.Figure 9
Efficiency comparison of (a) ACA and (b) bMVP.
(a)
(b)


### 3.3. Numerical Results of Deeply Optimized CG-FFT

To validate the accuracy and efficiency of the proposed method, the quasiperiodical structures of flat lens on substrate (εr = 2.2, μr = 1.0) with a copper ground (σ = 5.98 × 107 S/m) are considered, which consist of 15 × 15 and 33 × 33 Jerusalem crosses and generate 19561 and 85849 RWG unknowns, respectively, as shown in Figure 10. The flat lenses are illuminated by an orthogonal incident plane wave with polarization in the x direction, where the incident information areθ = π, φ = 0,E0x = 1, andE0y =E0z = 0.Figure 10
Geometry of flat lens on substrate and its meshing scheme.Since the MVP operation is the most time consuming part in CG-FFT, direct optimization from arithmetic programming point of view is a challenging topic and should be tested before applying to complicated structure. Therefore, both square and nonsquare matrices MVP are tested, from which a significant speedup ratio is observed comparing to the famous CUBLAS, as shown in Tables4 and 5. For this experiment on arithmetic realization, the GPU occupancy of CUBLAS and that of the proposed method is 25% and 75%, respectively, while the usage of shared memory is 288 B and 10240 B, respectively. This can be seen as a great enhancement in hardware performance, which makes a speedup ratio of 1.4 in average on single GPU.Table 4
Comparison in MVP for square matrix.
Matrix dimension
CUBLAS (ms)
This paper (ms)
4,000
1.58
1.12
5,000
2.50
1.77
6,000
3.43
2.54
7,000
4.63
3.46
8,000
5.65
4.07
9,000
7.17
5.19
10,000
8.83
6.48
11,000
10.75
7.94
12,000
13.28
9.53Table 5
Comparison in MVP for nonsquare matrix.
Matrix dimension
CUBLAS (ms)
This paper (ms)
1,700 × 9,000
1.69
1.10
1,800 × 9,000
1.58
1.17
1,900 × 9,000
1.71
1.22
2,900 × 15,000
4.28
2.95
3,000 × 15,000
4.00
2.98
3,100 × 15,000
4.69
3.16
4,900 × 25,000
11.48
8.24
5,000 × 25,000
11.19
8.44
5,200 × 25,000
11.29
8.72As is mentioned in Section2, the right-hand term in MoM of multilayered structure is different from radiating, transmitting, and scattering cases. With an orthogonal incident wave, the right-hand term can be calculated as below.(16)Vm=-jωε0∫TmE⃑incr⃑·f⃑mr⃑ds,E⃑incρ=E⃑0ejkxx+kyy+kzd+R~ejkxx+kyy-kzd,where the calculation of reflection coefficient refers to [37] andd in the exponential term is the thickness of substrate. The relative error of iteration in CG-FFT for the case with 15 × 15 elements is shown in Figure 11 and the current distribution in centre of the case with 33 × 33 elements is shown in Figure 12. To enhance the efficiency further, we have realized the fast calculation of multilayered Green’s functions on GPU/CPU heterogeneous platform [38], which will be implemented to the MPIE-MoM in future work.Figure 11
Relative error of iterations of the flat lens with 15 × 15 elements.Figure 12
Current distribution in centre of the flat lens with 33 × 33 elements.


## 3.1. Numerical Results of OpenMP-CUDA-MLFMM

As a benchmark result to verify the robustness of algorithms, the scattering from the PEC NASA Almond at 7 GHz for VV polarization and that from the PEC NASA Ogive at 9 GHz for HH polarization are considered. These two PEC surfaces are discretized into refined triangles with edge length of about 0.1λ, which produce 8416 and 15516 RWG unknowns, respectively. The monostatic RCS curves with the incident angle as argument are obtained by the proposed method, as shown in Figure 6. From the comparison with measurement results, very good agreement can be found. Note that the interpretations labeled as “mid VV,” “high VV,” and “VV FERM” refer to Figure 5 in [36], and those labeled as “HH,” “VV,” “HH Cicero,” and “VV Cicero” refer to Figure 9 in [36]. To demonstrate the efficiency performance of the proposed OpenMP-CUDA-MLFMM, we need to define some variables to record time consumption in different procedures, as shown in Table 3.Table 3
Variables definition of time consumption.
Variables
Description
t
M
Time for near matrix generation
t
S
Time for corresponding SpMVP
t
A
F
Time for aggregation at finest level
t
D
F
Time for disaggregation at finest level
t
A
S
Time for aggregation towards coarsest level
t
D
S
Time for disaggregation towards finest levelFigure 6
(a) The monostatic RCS at 7 GHz for VV polarization. (b) The monostatic RCS at 9 GHz for HH polarization.
(a)
(b)By sampling on different number of unknowns, an approximate curve of(tM;tS) and (tAF+tDF;tAS+tDS) can be fitted in a broad range, as shown in Figure 7. In can be clearly seen that the proposed method has higher efficiency than that reported in [16] based on the same Kepler GK110 architecture. It is worth pointing out that the translation is not included because both the proposed method and that in [16] have almost the same performance in this procedure.Figure 7
t
M and tS as function ofN with (a) single and (b) double precision; tAF+tDF and tAS+tDS as function ofN with (c) single and (d) double precision.
(a)
(b)
(c)
(d)


## 3.2. Numerical Results of GPU-Based ACA

To verify the robustness of the proposed GPU-based ACA, the bistatic RCS from a rotor with four blades and a depressed cylindrical object are considered, as shown in Figure8. Very good agreement can be seen between the proposed method and the traditional MoM on CPU. These two PEC surfaces are discretized into refined triangles with edge length of about 0.1λ, which produce 6636 and 4817 RWG unknowns, respectively.Figure 8
Bistatic RCS of (a) 4 blades rotor and (b) depressed cylindrical object.
(a)
(b)To demonstrate the efficiency performance of the proposed GPU-based ACA, the scattering from a PEC sphere with radius ofnλ (n= 1, 2, 3, and 4) is considered. The number of generated RWG unknowns is 4774, 18468, 40950, and 72120, respectively. In this paper, the threshold in formula (6) is selected as εACA=10-3. The efficiency comparison between GPU and CPU platform in both the bMVP operations and the whole ACA scheme with single and double precision are shown in Figure 9. It can be seen that the speedup ratio of matrix compression by using ACA can achieve about 50~100 and 25~50 for single and double precision, respectively. Meanwhile, the speedup ratio of bMVP can achieve about  10~30 and 6~17 for single and double precision, respectively.Figure 9
Efficiency comparison of (a) ACA and (b) bMVP.
(a)
(b)


## 3.3. Numerical Results of Deeply Optimized CG-FFT

To validate the accuracy and efficiency of the proposed method, the quasiperiodical structures of flat lens on substrate (εr = 2.2, μr = 1.0) with a copper ground (σ = 5.98 × 107 S/m) are considered, which consist of 15 × 15 and 33 × 33 Jerusalem crosses and generate 19561 and 85849 RWG unknowns, respectively, as shown in Figure 10. The flat lenses are illuminated by an orthogonal incident plane wave with polarization in the x direction, where the incident information areθ = π, φ = 0,E0x = 1, andE0y =E0z = 0.Figure 10
Geometry of flat lens on substrate and its meshing scheme.Since the MVP operation is the most time consuming part in CG-FFT, direct optimization from arithmetic programming point of view is a challenging topic and should be tested before applying to complicated structure. Therefore, both square and nonsquare matrices MVP are tested, from which a significant speedup ratio is observed comparing to the famous CUBLAS, as shown in Tables4 and 5. For this experiment on arithmetic realization, the GPU occupancy of CUBLAS and that of the proposed method is 25% and 75%, respectively, while the usage of shared memory is 288 B and 10240 B, respectively. This can be seen as a great enhancement in hardware performance, which makes a speedup ratio of 1.4 in average on single GPU.Table 4
Comparison in MVP for square matrix.
Matrix dimension
CUBLAS (ms)
This paper (ms)
4,000
1.58
1.12
5,000
2.50
1.77
6,000
3.43
2.54
7,000
4.63
3.46
8,000
5.65
4.07
9,000
7.17
5.19
10,000
8.83
6.48
11,000
10.75
7.94
12,000
13.28
9.53Table 5
Comparison in MVP for nonsquare matrix.
Matrix dimension
CUBLAS (ms)
This paper (ms)
1,700 × 9,000
1.69
1.10
1,800 × 9,000
1.58
1.17
1,900 × 9,000
1.71
1.22
2,900 × 15,000
4.28
2.95
3,000 × 15,000
4.00
2.98
3,100 × 15,000
4.69
3.16
4,900 × 25,000
11.48
8.24
5,000 × 25,000
11.19
8.44
5,200 × 25,000
11.29
8.72As is mentioned in Section2, the right-hand term in MoM of multilayered structure is different from radiating, transmitting, and scattering cases. With an orthogonal incident wave, the right-hand term can be calculated as below.(16)Vm=-jωε0∫TmE⃑incr⃑·f⃑mr⃑ds,E⃑incρ=E⃑0ejkxx+kyy+kzd+R~ejkxx+kyy-kzd,where the calculation of reflection coefficient refers to [37] andd in the exponential term is the thickness of substrate. The relative error of iteration in CG-FFT for the case with 15 × 15 elements is shown in Figure 11 and the current distribution in centre of the case with 33 × 33 elements is shown in Figure 12. To enhance the efficiency further, we have realized the fast calculation of multilayered Green’s functions on GPU/CPU heterogeneous platform [38], which will be implemented to the MPIE-MoM in future work.Figure 11
Relative error of iterations of the flat lens with 15 × 15 elements.Figure 12
Current distribution in centre of the flat lens with 33 × 33 elements.


## 4. Conclusions

In this paper, our recent works on high performance computing based on GPU/CPU heterogeneous platform and its application to EM scattering problems and planar multilayered medium structures are introduced. There are three typical applications analyzed, including a novel implementation of OpenMP-CUDA-MLFMM, a developed GPU-based ACA method, and a deeply optimized CG-FFT method. Comparing to the MLFMM, the ACA method starts to accelerate computing with Green’s functions; however, it is still worth mentioning that the data storage strategy of ACA needs to be further developed. This paper is to explain the coordination between GPU and CPU on a heterogeneous platform. In addition, the methodology of optimization through hardware instructions is also adopted. Following the trends of both parallel and distributed computing, the OpenMP-CUDA-MLFMM and ACA are realized on a 2-GPU platform. With fruitful numerical examples for the three kinds of algorithms, we can find obvious enhancement in efficiency. As the development in hardware world, it is worthy keeping an eye on different kinds of modules, including GPU, FPGA, and DSP. Meanwhile, it is also worthy keeping one on investigating and understanding the operating mechanism of heterogeneous computing.



---

*Source: 9173062-2017-07-30.xml*