Bogaert, Ignace

De Zutter, Daniël

Fostier, Jan

Michiels, Bart

English

Ghent University Academic Bibliography

Weak Scaling Analysis for the Parallel MLFMA:The Quest for Ever Larger SimulationsBart Michiels1, Jan Fostier1, Ignace Bogaert1, Danie¨l De Zutter11Department of Information Technology (INTEC), Ghent University, BelgiumAbstract—In the past decade, several data distribution strate-gies for the parallel, distributed-memory Multilevel Fast Multi-pole Algorithm (MLFMA) have been introduced. The commongoal is to distribute the computations in the MLFMA uniformlyover the different parallel processes, while minimizing dependen-cies and data communication between them. However, as clusterswith thousands of CPU cores are becoming increasingly available,the asymptotic behavior of algorithms for a very large number ofprocesses P and problem size N (weak scaling analysis) starts toplay a dominant role. Parallel MLFMA implementations basedon a hierarchical distribution of boxes and radiation patternsampling points that exhibit excellent weak scaling behaviorappear to be attractive candidates to tackle even larger problemson future hardware clusters.Index Terms—parallel MLFMA, weak scaling analysisI. INTRODUCTIONThe Multilevel Fast Multipole Algorithm (MLFMA) is a fastmatrix-vector multiplication scheme, used during the iterativesolution of a Method of Moments (MoM) discretization oflarge boundary integral equation problems. Its computationalcomplexity is O(N logN) [1], with N the number of un-knowns. State-of-the-art sequential MLFMA implementationscan handle over a million of unknowns on a single worksta-tion. However, several real-life, electrically large geometriesrequire discretization into hundreds of millions, if not billionsof unknowns. In order to handle those effectively, parallel,distributed-memory implementations are developed in orderto take advantage of the aggregated memory and computecapacity of clusters.In [2], the weak scaling behavior of a number of existingpartitioning schemes was examined. Weak scaling impliesthat larger problems can be solved on proportionally largercomputational clusters without loss of efficiency. In a weakscaling setup, the number of processes P scales linearly withN , i.e., P = O(N) and hence, the number of unknowns perprocess is kept fixed. Given a sequential complexity for theMLFMA of O(N logN), this implies that the complexity pernode should not exceed O(1) per level, or hence, O(logN)in total.None of the existing schemes for the parallel MLFMA(spatial partitioning (SP), hybrid partitioning (HyP) [3] andhierarchical partitioning (S-HiP) [4], [5]) were found to exhibitgood weak scaling behavior. In all cases, the worst-casecomputational complexity per node and per level exceeds O(1)(see entries 1-3 of Table I). As a consequence, when consid-ering larger problems on larger clusters, the load imbalanceTABLE IOVERVIEW OF THE DIFFERENT PARTITIONING STRATEGIES AND THEWORST-CASE COMPUTATIONAL COMPLEXITY FOR A NODE PER LEVEL.partitioning worst-case complexitystrategy per node and per levelspatial (SP) O(N)hybrid (HyP) O(√N)hierarchical - strip (S-HiP) O(√N)hierarchical - block (B-HiP) O(1)between processes will grow, leaving many processes idlewhile only a few of them are performing useful calculations.The major difficulty in the parallelization of the MLFMAis the fact that the data structures and the relative contributionfrom different types of computations differ from level to levelin the MLFMA tree. At the top levels, there are only fewboxes, each box containing a large number of radiation patternsampling points, whereas at the lower levels, there are manyboxes, each box containing much smaller radiation patterns.Existing data distribution schemes (SP, HyP and S-HiP) failto distribute the load evenly on all levels (for sufficiently highN and P ). For the SP, the load imbalance becomes alreadyapparent for a very modest number of parallel processes (e.g.at 16 processes). For the more advanced S-HiP scheme, theload imbalance only manifests itself at a very high number ofparallel processes and very large problem sizes.We proposed a modification of the hierarchical scheme,in which the sampling points on the sphere are uniformlysampled in elevation θ and azimuth φ [6], and distributed inboth θ and φ [2]. This scheme, denoted by B-HiP (see lastentry of Table I), was found to exhibit excellent weak scalingbehavior (i.e., worst-case complexity of O(1) per node andper level [7]), at the cost of being significantly more complexto implement.Nowadays, computational clusters with thousands of CPUare available, enabling the solution of problems with hundredsof millions or even billions of unknowns. For such extremelylarge-scale problems and high number of parallel processes,the load imbalance in the SP, HyP and even the S-HiPscheme are becoming apparent. In the next section, we analyzesuch large-scale problem and indicate where the bottlenecksemerge. At the time of conference, the accurate solution ofa problem involving billions of unknowns using thousands ofCPU cores will be presented.II. NUMERICAL EXAMPLEWe examine the behavior of all three previously existingpartitioning strategies (SP, HyP and S-HiP) for a large, butnowadays realistic problem size N and a high, but againrealistic number of parallel processes P . For each of theschemes, the load imbalance that causes bad performance ispinpointed. We emphasize that these bottlenecks will becomeeven more pronounced when N and P are increased evenfurther, as can be expected with future hardware and problems.Consider a perfectly electrically conducting (PEC) cubegeometry with an edge size of 1638λ. Using a λ/10 dis-cretization, this problem contains slightly over 1.2 billionof unknowns. Table II shows, for the different levels of theMLFMA tree, the sampling rate of the radiation patterns andthe number of boxes at the level of tree. Assume that wewish to solve this problem using 4096 parallel processes. Thisresults in roughly 300 000 unknowns per process.If this problem were to be solved using the spatial partition-ing (SP) technique, then the boxes on all levels would haveto be partitioned evenly among all parallel processes. Clearly,this is only possible for the lowest levels. For levels 6 andhigher, the number of boxes is even lower than the numberof parallel processes. At these levels, the load balancing willbe particularly problematic as there are fewer boxes thanprocesses, leaving many processes idling. At level 9, only 56out of 4096 processes would be active. However, the amountof work at this level is roughly the same as on all other levels.This explains why the SP achieves only very limited efficiency.The same reasoning can be applied to the hybrid (HyP) schemeat the transition level (i.e., the highest level that is still usingspatial partitioning, typically the middle level).If the hierarchical partitioning technique with a one-dimensional partitioning of radiation patterns (S-HiP) were tobe used (e.g. only in the θ or φ direction), then at level 9,the 2871 sampling directions along θ or the 5744 samplingdirections along φ would have to be distributed uniformlyamong the 4096 processes. Again, it is impossible to obtain auniform partitioning, however, the amount of load imbalance isclearly less severe than in the SP or HyP case. Nevertheless,the number of sampling points in either θ or φ direction ata top level only grows as O(√N) whereas the number ofparallel processes grows as P = O(N). This means that theload imbalance will effectively become more severe for higherN and P .This analysis also explains why smaller simulations (tensor hundreds of millions of unknowns) using a (for currentstandards) moderate number of parallel processes (e.g. 128)can effectively be performed with high efficiency using theHyP or S-HiP approach. Whereas the breakdown for the SPtechnique becomes apparent already at 16 parallel processes,the breakdown for the S-HiP scheme appears only when thenumber of processes exceeds the number of sampling pointsin the θ or φ direction. This occurs only when consideringthe very largest problems with over thousands of parallelprocesses.TABLE IISAMPLING RATE AND NUMBER OF BOXES PER LEVEL OF THE MLFMATREE FOR A CUBE GEOMETRY WITH A SIZE OF 1638 λ.MLFMA box size radiation pattern number oflevel sampling rate (θ × φ) boxes0 0.5λ 17× 36 1 6098 2661 1λ 26× 54 4 024 5682 2λ 40× 82 1 003 6883 4λ 66× 134 249 6984 8λ 115× 232 62 4265 16λ 209× 420 15 6086 32λ 392× 786 3 7527 64λ 751× 1504 8668 128λ 1462× 2924 2189 256λ 2873× 5748 5610 512λ 5683× 11368 811 1024λ 11286× 22574 1III. CONCLUSIONWe reviewed several schemes (SP, HyP, S-HiP, B-HiP) forthe distributed-memory parallel MLFMA, and reviewed theirweak scaling behaviour by asymptotic analysis and througha numerical example for finite N and P . Only the B-HiPscheme exhibits good weak scaling behaviour. The B-HiPscheme attains a per node complexity of O(logN) when usingP = O(N) parallel processes. As future clusters will likelyincorporate more and more parallelism, the B-HiP schemeappears to be the most attractive approach to handling evenlarger problems.ACKNOWLEDGMENTThe computational resources (Stevin Supercomputer Infras-tructure) and services used in this work were provided byGhent University, the Hercules Foundation and the FlemishGovernment department EWI. The work of B. Michiels wassupported by a doctoral grant from the Special ResearchFund (BOF) at Ghent University. The work of I. Bogaert wassupported by a post-doctoral grant from Research Foundation-Flanders (FWO-Vlaanderen).REFERENCES[1] W.C. Chew, J. Lin, E. Michielssen and J. Song, Fast and Efficient algo-rithms in Computational Electromagnetics. Boston: Artech House, 2001.[2] B. Michiels, J. Fostier, I. Bogaert, and D. De Zutter, “Weak ScalabilityAnalysis of the Distributed-Memory MLFMA,” accepted for publicationin IEEE Transactions on Antennas and Propagation, 2013.[3] S. Velamparambil and W.C. Chew, “Analysis and Performance of aDistributed Memory Multilevel Fast Multipole Algorithm”, IEEE Trans-actions on Antennas and Propagation, vol. 53, no. 8, pp. 2719–2727,2005.[4] O¨. Ergu¨l and L. Gu¨rel, “A Hierarchical Partitioning Strategy for anEfficient Parallelization of the MultiLevel Fast Multipole Algorithm”,IEEE Tra nsactions on Antennas and Propagation, vol. 57, no. 6, pp.1740–1750, June 2009.[5] O¨. Ergu¨l, and L. Gu¨rel, “Hierarchical Parallelisation Strategy for Mul-tilevel Fast Multipole Algorithm in Computational Electromagnetics,”Electronics Letters, vol. 44, no. 1, pp. 3–4, Jan. 2008.[6] J. Sarvas, “Performing Interpolation and Anterpolation entirely by FastFourier Transform in the 3-D Multilevel Fast Multipole Algorithm”,SIAM Journal on Numerical Analysis, 41(6):2180–2196, 2003.[7] J. Fostier, and F. Olyslager, “Provably Scalable Parallel Multilevel FastMultipole Algorithm,” Electronics Letters, vol. 44, no. 19, pp. 1111–1112,Sept. 2008.

Weak scaling analysis for the parallel MLFMA: the quest for ever larger simulations

https://biblio.ugent.be/publication/5840444/file/5840448.pdf

Weak scaling analysis for the parallel MLFMA: the quest for ever larger simulations

Abstract

Similar works

Full text

Available Versions

Ghent University Academic Bibliography