2SWLPL]DWLRQRI6HULDODQG3DUDOOHO&RPPXQLFDWLRQV
IRU3DUDOOHO*HRPHWULF0XOWLJULG0HWKRG୰ ᓥ ◊ ࿃ᮾி኱Ꮫ᝟ሗᇶ┙ࢭࣥࢱ࣮
ᮏ✏䛷䛿䠈 ᖺ  ᭶䠈 ᖺ  ᭶䛻ᐇ᪋䛥䜜䛯኱つᶍ +3& 䝏䝱䝺䞁䝆䛾⤖ᯝ䜢ሗ࿌䛩䜛䚹ᮏ኱つᶍ +3& 䝏䝱䝺䞁䝆䛾ᡂᯝ䜢䜎䛸䜑䛯ㄽᩥ䛂.HQJR 1DNDMLPD 2SWLPL]DWLRQ RI 6HULDO DQG 3DUDOOHO&RPPXQLFDWLRQVIRU3DUDOOHO*HRPHWULF0XOWLJULG0HWKRG䛃䛿䠈 ᖺ  ᭶䛻᪂➉䠄ྎ‴䠅䛷㛤ദ䛥䜜䛯 ,(((,&3$'6 䠄WK,QWHUQDWLRQDO&RQIHUHQFHIRU3DUDOOHODQG'LVWULEXWHG6\VWHPV䠅䛻᥇ᢥ䛥
䜜䠈3URFHHGLQJV㸦,((( 'LJLWDO /LEUDU\㸧࡟ᥖ㍕ࡉࢀࡿ࡜࡜ࡶ࡟㸪%HVW 3DSHU $ZDUG ࢆཷ㈹ࡋࡓࠋ
,&3$'6 ࡣྜィ  ௳ࡢᢞ✏ࡀ࠶ࡾ㸪ࡑࡢ࠺ࡕ  ௳ࡀ᥇ᢥࡉࢀ㸦᥇ᢥ⋡ 㸧᭦࡟㸪ࡑࡢ࠺ࡕ  ௳ࡀ %HVW3DSHU&DQGLGDWHV ࡜ࡋ࡚㑅ฟࡉࢀ㸪᭱⤊ⓗ࡟ᮏㄽᩥࡀ %HVW3DSHU$ZDUG ࢆཷ㈹ࡋࡓࠋᮏ✏ࡣฟ∧ඖ࡛࠶ࡿ ,((( ࡢチྍࢆཷࡅḟ࣮࣌ࢪ௨㝆࡟ୖグㄽᩥࡢ඲ᩥࢆᥖ㍕ࡍࡿࡶࡢ࡛࠶ࡿࠋ
5HSULQWHG ZLWK SHUPLVVLRQ IURP ,((( IXOO FUHGLW WR WKH RULJLQDO VRXUFH .HQJR 1DNDMLPD2SWLPL]DWLRQ RI 6HULDO DQG 3DUDOOHO &RPPXQLFDWLRQV IRU 3DUDOOHO *HRPHWULF 0XOWLJULG 0HWKRG
3URFHHGLQJV RI WKH WK ,((( ,QWHUQDWLRQDO &RQIHUHQFH IRU 3DUDOOHO DQG 'LVWULEXWHG 6\VWHPV
,&3$'6  IROORZHGE\WKH,(((FRS\ULJKWOLQHF>@,((( http://www.icpads.org/

スーパーコンピューティングニュース- 14 -

Vol.17,No.6

2015

Optimization of Serial and Parallel Communications
for Parallel Geometric Multigrid Method
Kengo Nakajima
Information Technology Center, The University of Tokyo
2-11-16 Yayoi, Bunkyo-ku, Tokyo 113-8658, Japan
nakajima@cc.u-tokyo.ac.jp
Abstract— The parallel multigrid method is expected to play
an important role in large-scale scientific computing on postpeta/exa-scale supercomputer systems, and it also includes serial
and parallel communication processes which are generally
expensive. In the present work, new format for sparse matrix
storage based on sliced Ellpack-Itpack (ELL) format is proposedfor optimization of serial communication in data transfer
through memories, and hierarchical coarse grid aggregation(hCGA) is introduced for optimization of parallel
communication by message passing. The proposed methods are
implemented for pGW3D-FVM, a parallel code for 3D
groundwater flow simulations using the multigrid method, andthe robustness and performance of the code was evaluated on up
to 4,096 nodes (65,536 cores) of the Fujistu FX10 supercomputer
system at the University of Tokyo. The parallel multigrid solverusing the sliced ELL format provided performance improvement
in both weak scaling (25%–31%) and strong scaling (9%–22%)
compared to the code using the original ELL format. Moreover,
hCGA provided excellent performance improvement in both
weak scaling (1.61 times) and strong scaling (6.27 times) for flat
MPI parallel programming model.
Keywords—parallel computing; iterative solvers; multigrid;communication; matrix storage format; multicore

I. INTRODUCTION
A multigrid is a scalable method for solving linear
equations and preconditioning Krylov iterative linear solvers,
and is especially suitable for large-scale problems. The
parallel multigrid method is expected to be one of thepowerful tools on post-peta/exa-scale systems. Recently,
HPCG (High Performance Conjugate Gradients) [1] wasproposed as a new benchmark for evaluation of the practical
performance of supercomputer systems. HPCG solves sparse
matrices derived from finite-element application using
conjugate gradient linear solver (CG) preconditioned withmultigrid method.
In previous works by the author [2,3,4], OpenMP/MPI
hybrid parallel programming models were implemented for
pGW3D-FVM, a 3D finite-volume simulation code forgroundwater flow problems through heterogeneous porous
media, by using parallel conjugate gradient (CG) solver withmultigrid preconditioner (MGCG). The performance and therobustness of the developed code were evaluated on multicore
clusters, such as the T2K Open Supercomputer (T2K/Tokyo)
and the Fujitsu FX10 System (Oakleaf-FX) at the Universityof Tokyo [5], by using up to 4,096 nodes (65,536 cores) for
both weak and strong scaling computations.

スーパーコンピューティングニュースIt is well-known that convergence of the solver at thecoarsest level of the multigrid cycle (coarse grid solver)
strongly affects convergence of the entire process of multigrid[2,6,7,8]. The coarse grid aggregation (CGA) proposed in [3]improves the performance and the robustness of multigrid
procedures with large numbers of MPI processes. In [4], theeffect of a format of sparse matrix storage on the performanceof MGCG was evaluated. The Ellpack-Itpack (ELL) format
was applied to pGW3D-FVM, and it provided excellent
improvement of memory access throughput, and the MGCG
solver using the ELL format with CGA showed excellent
scalable performance and robustness. The performance
improvement from the original solver with the compressedrow storage (CRS) format in [3] to the new one with ELLCGA at 4,096 nodes of the Fujitsu FX10 System (OakleafFX) was 13%–35% for weak scaling, and 40%–70% forstrong scaling.
The parallel multigrid method and MGCG include both ofserial and parallel communication processes which are
generally expensive. The serial communication is the datatransfers through memory hierarchies of each processor.
Because MGCG solver with sparse coefficient matrices is amemory-bound process, serial communication is a serious
problem. In the present work, further modification is
introduced into the ELL format of the coefficient matrices.
The parallel communication is by message passing between
computing nodes through the network using MPI. Hierarchicalcoarse grid aggregation (hCGA) is introduced for optimizationof the parallel communication in the present work.
The rest of this paper is organized as follows. In Section II,
an overview of the target hardware and application is provided,and results of the previous works are briefly overviewed. In
Section III, we give a summary of the method for optimizationof serial and parallel communications. Finally, Section IV
provides the results of computations, and final remarks areoffered in Sections V.
II. HARDWARE ENVIRONMENT AND TARGET APPLICATION
A. Hardware Environment
The Fujitsu FX10 system at the University of Tokyo(Oakleaf-FX) [9] is Fujitsu’s PRIMEHPC FX10 massively
parallel supercomputer with a peak performance of 1.13
PFLOPS. The Oakleaf-FX system consists of 4,800
computing nodes of SPARC64™ IXfx with 16 cores (1.848
GHz) [9]. The entire system consists of 76,800 cores and 154- 15 -

Vol.17,No.6

2015

TB memory. Each core has a 64 KB L1 instruction/data cache.
A 12 MB L2 cache is shared by 16 cores on each node. On the
SPARC64™ IXfx, each of the 16 cores accesses memory in a
uniform manner. Nodes are connected via a 6-dimensional
mesh/torus interconnect, called Tofu. Although users can
specify the topology of the network on Fujitsu FX10, thiscapability was not used in the present work.
B. pGW3D-FVM with MGCG Solver
pGW3D-FVM, a parallel simulation code based on thefinite-volume method (FVM), solves groundwater flow
problems through saturated heterogeneous porous media (Fig.
1). The problem is described by the following Poisson
equation and boundary condition (1):  O x, y , z I  q,I  0at z  
z max        (1)where Idenotes the potential of the water head, O(x,y,z)
describes the water conductivity, and q is the value of thevolumetric flux of each finite-volume mesh and is set to a
uniform value (=1.0) in this work.
(a)domain. The reordering of meshes in each domain allows theconstruction of local operations without global dependency
for achieving parallel IC operations in multigrid processes.
Cuthill-McKee (CMK), Reverse Cuthill-McKee (RCM), andRCM with cyclic multicoloring (CM-RCM) [2,3,4] are
implemented.
C. CGAIn [2], multigrid processes by V-cycle were applied asshown in Fig. 2. Multigrid operations at each level are done ina parallel manner, but the operations at the coarsest level(coarse grid solver) are executed on a single core by gathering
the information of entire processes. The total number ofmeshes at the coarsest level is equal to the number of MPI
processes. Moreover, overhead of parallel communication issignificant at coarser level of the V-cycle.
Fine
Level=1
Level=2(b)

Level=m-3
Level=m-2
Level=m-1
Level=mMesh # foreach MPI= 1

Fig.1 Example of groundwater flow through heterogeneous porous media. (a)Distribution of water conductivity; (b) Streamlines

A heterogeneous distribution of water conductivity in eachmesh is calculated by a sequential Gauss algorithm, which is
widely used in the area of geostatistics [10]. The minimum
and maximum values of water conductivity are 10-5 and 105,
respectively, with an average value of 1.0. Each mesh is acube, and distribution of the meshes is structured as finitedifference-type voxels. In the present work, an entire model
consists of clusters of small models with 1283 meshes withsame pattern [2,3,4]. The conjugate gradient (CG) solver withmultigrid preconditioner (MGCG) [2,3,4] was applied forsolving Poisson’s equations with symmetric positive definite(SPD) coefficient matrices derived by pGW3D-FVM. A very
simple geometric multigrid with a V-cycle algorithm is applied,
where 8 children form 1 parent mesh in an isotropic mannerfor structured finite-difference-type voxels. The level of thefinest grid is set to 1, and the level is numbered from the finest
to the coarsest grid, at which the number of meshes is 1 ateach MPI process. Incomplete Cholesky factorization without
fill-ins (IC(0)) is adopted as a smoothing operator of multigrid
process for ill-conditioned problems. The additive Schwarz
domain decomposition (ASDD) for overlapped regions [11] isintroduced for stabilization of the block-Jacobi-type localized
procedure of parallel IC(0). The pGW3D-FVM code isparallelized by domain decomposition using MPI [2,3,4]. In
the OpenMP/MPI hybrid parallel programming model,
multithreading by OpenMP is applied to each partitioned

スーパーコンピューティングニュース

CoarseCoarse grid solver on a single
core (further multigrid)Fig.2. Original parallel multigrid method with V-cycle (restriction) [2]Fine
Level=1
Level=2

Level=m-3
Level=m-2

Coarse grid solver on asingle MPI process (multithreaded, further
multigrid)

CoarseFig.3. Procedures of coarse grid aggregation (CGA), where information ofeach MPI process is gathered in a single MPI process for computation atlevel=m-2

In [3], CGA was proposed, where operations for
aggregation/disaggregation of MPI processes in Fig. 2 aredone at a finer level. If we switch to a coarse grid solver at afiner level, more robust convergence and reduction ofcommunication overhead are expected, even though the sizeof the coarse grid problem is larger than that of the original
configuration. Furthermore, the coarse grid solver is multi-- 16 -

Vol.17,No.6

2015

ª1
«1
«
«4
«
«0
«¬1

3 0 0 0º
2 5 0 0»»1 3 0 0»»3 7 4 0»0 0 0 5»¼

(a)

(b)

1

3

1

3

0

1

2

5

1

2

5

4

1

3

4

1

3

3

7

4

3

7

4

1

5

1

5

0

Fig. 4. Formats of sparse matrix storage. (a) Compressed row storage (CRS);(b) Ellpack-Itpack (ELL)

E. Results in the Previous Works
The performance of the developed code was evaluated on 8
to 4,096 nodes of the Fujitsu FX10. The following three typesof OpenMP/MPI hybrid parallel programming tablemodels
were applied, and the results were compared with those of flat
MPI:
x
x
x

Hybrid 4u4 (HB 4u4): Four OpenMP threads for each
MPI process, four MPI processes on each nodeHybrid 8u2 (HB 8u2): Eight OpenMP threads for each
MPI process, two MPI processes on each nodeHybrid 16u1 (HB 16u1): Sixteen OpenMP threads foreach MPI process, a single MPI process on each node.

The number of finite-volume meshes per core was 262,144(=643). The performance of solvers using CRS, ELL, and ELL
with CGA (ELL-CGA) were evaluated. Figure 5 compares
the performances of HB 8u2 at 4,096 nodes (65,536 cores). A
2-digit number in each label after a colon on the transverse
axis of Fig.5 indicates the number of iterations untilconvergence. Because the coarse grid solver was performed
on a single MPI process in ELL-CGA, the convergence
provided by IC(0) smoothing is much more robust than that by
localized IC(0) smoothing [3].

スーパーコンピューティングニュースRest
Coarse Grid Solver
MPI_Allgather
MPI_Isend/Irecv/Allreduce

sec.

10.0

5.0

0.0

ELL-CGA, ELL-CGA, ELL-CGA, ELL: 65, CRS: 66,
lev=6: 51 lev=7: 55 lev=8: 60 (NO CGA) (NO CGA)

Fig. 5. Performance of MGCG solver on Fujitsu FX10 using 4,096 nodes(65,536 cores), total problem size: 17,179,869,184 meshes, comparison of
CRS and ELL, effect of switching level for coarse grid solver in ELL-CGA,
HB 8u215.0

Flat MPI:ELL
HB 4x4:ELL-CGA
HB 8x2:ELL-CGA
HB 16x1:ELL-CGA

12.5

sec.

D. CRS and ELL
Generally, computations with sparse matrices, such as
multigrid and MGCG, are memory-bound processes, becauseof the indirect memory accesses. Various types of storage
formats have been proposed. The compressed row storage(CRS) format is the most popular and widely used because of
its flexibility. It stores only non-zero components of sparse
matrices, as shown in Fig. 4(a). In the Ellpack-Itpack (ELL)
format, the number of non-zero components of each row is set
to that of the longest non-zero entry row of the matrix, as
shown in Fig. 4(b). This format allows one to achieve better
performance for memory access than CRS, but introduces
extra computations and memory requirements, since somerows are zero-padded, as shown in Fig. 4(b).

15.0

10.0

7.5

5.0

100

1000

10000

100000

CORE#

Fig. 6. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), weak scaling: 262,144 (=643) meshes/core, max. total
problem size: 17,179,869,184 meshes, elapsed computation time for MGCG,
switching to coarse grid solver at level=7 for ELL-CGA cases
4.00

Flat MPI
HB 8x2

HB 4x4
HB 16x1

3.00

sec.

threaded by OpenMP and uses all cores on each MPI process,although only a single core on each node was utilized in theoriginal method [2]. Figure 3 shows procedures of CGA,
where information of each MPI process is gathered in a single
MPI process for computation at level=m-2. Thus, the stage ofthe coarse grid solver starts earlier than that starts in theoriginal case. CGA is not applied to flat MPI.

2.00

1.00

0.00

1024 2048 4096 8192 16384 32768 49152 65536
CORE#

Fig. 7. Computation time for the coarse grid solver of MGCG solver on
Fujitsu FX10 using 4,096 nodes (65,536 cores), total problem size:
17,179,869,184 meshes for ELL-CGA cases

The optimum parameter is “level=7” in this case, which has
been determined through empirical parameter studies.
Improvement of the performance from CRS to ELL-CGA is
13%–35% at 4,096 nodes [3]. Figure 6 provides the results of
weak scaling. Switching at “level=7” was applied to all ELLCGA cases. Although each of HB 4u4, HB 8u2 and HB 16u1
provides excellent scalability, HB 8u2 with ELL-CGA- 17 -

Vol.17,No.6

2015

provides the best performance at 4,096 nodes. Flat MPI withmore than 104 cores cannot achieve scalability.
Figure 7 shows that the computation time spent for thecoarse grid solver is getting larger for flat MPI with morethan 104 cores. This is because the initial problem size for thecoarse grid solver is larger, and coarse grid problems are
solved by a single core in flat MPI cases. In contrast, the sizeof the coarse grid problem for HB 16u1 is smaller and the
problem is solved by 16 threads in HB 16u1. Computationalamount of coarse grid solver for each core of flat MPI is 256(=16u16) times as large as that of HB 16u1.
III. OPTIMIZATION OF SERIAL AND PARALLELCOMMUNICATIONS

Fig.9. Procedures for backward substitution by original ELL format [4],Cache is not well-utilized for loops of pure internal meshes

F. Sliced ELL for Serial Communication
In the previous work [4], ELL was introduced as a formatof sparse matrix storage. Although structured meshes are usedin pGW3D-FVM, utilization of ELL is not straightforward. In
pGW3D-FVM, each of the diagonal, lower triangular, andupper triangular components of the coefficient matrices are
separately stored in different arrays for IC(0) operations.
(a)External
MeshesInternal Meshes
on Boundary
Pure Internal
Meshesdo icol= NHYP(lev), 1, -1
if (mod(icol,2).eq.1) then!$omp parallel do private (ip,icel,j,SW) for Boundary Cells
do ip= 1, PEsmpTOT
do icel= SMPindex(icol-1,ip,lev)+1, SMPindex(icol,ip,lev)
SW= 0.0d0
do j= 1, 6
SW= SW + AUnew(j,icel)*Rmg(IAUnew(j,icel))
enddo
Rmg(icel)= Rmg(icel) - SW*DDmg(icel)
enddo
enddo
else
!$omp parallel do private (ip,icel,j,SW) for Pure Internal Cells
do ip= 1, PEsmpTOT
do icel= SMPindex(icol-1,ip,lev)+1, SMPindex(icol,ip,lev)
SW= 0.0d0
do j= 1, 3
SW= SW + AUnew(j,icel)*Rmg(IAUnew(j,icel))
enddo
Rmg(icel)= Rmg(icel) - SW*DDmg(icel)
enddo
enddo
endif
IAUnew (6,N), AUnew (6,N)
enddo

do icol= NHYP(lev), 1, -1
if (mod(icol,2).eq.1) then!$omp parallel do private (ip,icel,j,SW) for Boundary Cells
do ip= 1, PEsmpTOT
do icel= SMPindex(icol-1,ip,lev)+1, SMPindex(icol,ip,lev)
SW= 0.0d0
do j= 1, 6
SW= SW + AUnew6(j,icel)*Rmg(IAUnew6(j,icel))
enddo
Rmg(icel)= Rmg(icel) - SW*DDmg(icel)
enddo
enddo
else
!$omp parallel do private (ip,icel,j,SW) for Pure Internal Cells
do ip= 1, PEsmpTOT
do icel= SMPindex(icol-1,ip,lev)+1, SMPindex(icol,ip,lev)
SW= 0.0d0
do j= 1, 3
SW= SW + AUnew3(j,icel)*Rmg(IAUnew3(j,icel))
enddo
Rmg(icel)= Rmg(icel) - SW*DDmg(icel)
enddo
IAUnew3(3,N), AUnew3(3,N)
enddo
endif
IAUnew6(6,N), AUnew6(6,N)
enddo

Fig.10. Procedures for backward substitution by sliced ELL format,
IAUnew/AUnew are stored in separate arrays for boundary meshes and pureinternal meshes, Cache is well-utilized for loops of pure internal meshes

Figure 9 shows the procedures of backward substitution ofIC(0) smoothing for ELL-CGA [3]. Arrays for column ID’s
(IAUnew(6,N)), and coefficients (AUnew(6,N)) for upper
triangular components of IC(0) matrix are used for both of
䖃 Internal(b)pure internal meshes and internal meshes on boundary(lower)(boundary meshes) in Fig.8(b). At each level of Cuthill䖃 Internal
McKee (CMK) reordering (icol in Fig.9), pure internal(upper)meshes and boundary meshes are calculated separately.
z
Because reordering on a same level does not affecty䖃 Externalconvergence [2,3,4], computations related to upper triangular(upper)components such as the backward substitution process inInternal Meshes
Pure InternalIC(0) smoothing for boundary meshes (with up to 6 upperx
on Boundary
MeshesFig. 8. Local data structure of pGW3D-FVM. (a) Internal/external meshes; (b)components) are done first, then pure internal meshes (with up
Upper/lower components of internal meshesto 3) are calculated. This implementation of ELL provides
25% improvement of performance compared to CRS for HB
In Fig. 8, (a) and (b) provide the local data structure of 8u2 at 4,096 nodes, as shown in Fig.5. But, this approach is
pGW3D-FVM. The numbering of external meshes on each not necessarily efficient. IAUnew/AUnew(4,N)-(6,N)
distributed mesh starts after all internal meshes in Fig. 8(a) are not used for calculation of pure internal meshes, although
are numbered. In the structured mesh in pGW3D-FVM, the these are stored on cache. In the present work, improved
initial numbering is lexicographical. Therefore, each pure version of ELL (Sliced ELL) is introduced, whereinternal mesh has three lower and three upper triangular IAUnew/AUnew for pure internal meshes and boundary
components. But, an internal mesh on a domain boundary meshes
arestoredinseparatearrays,may have up to six upper triangular components, as shown in
IAUnew3/AUnew3(3,N) and IAUnew6/AUnew6(6,N),
Fig. 8(b). In [4], the ELL format is applied to pGW3D-FVM
as shown in Fig.10. This idea is very similar to Sliced
with Cuthill-McKee (CMK) reordering, because CMK does
ELLPACK proposed in [12]. In the Sliced ELLPACK, a sparsenot change the inequality relationship of the ID of each meshmatrix is first divided into sub-matrices (slices) consisting of S
for the structured meshes in pGW3D-FVM. In the ELL format,
non-zero off-diagonal components (S=1,…, N), and each slice
2D-type arrays for matrices (e.g. AMAT(6,N)) are applied.
is then stored in the ELL format. Thus, the redundancy

スーパーコンピューティングニュース- 18 -

Vol.17,No.6

2015

inherent in the ELL format is eliminated, and cache is wellutilized in loops for pure internal meshes. In the present work,
we have two slices for boundary meshes and pure internal
meshes, as shown in Fig.11. CMK reordering is also applied
for the new format of matrix storage based on sliced ELL. The
effect of formats of sparse matrix storage was evaluated using
a single node (16 cores) of Fujitsu FX10 for flat MPI. The
number of meshes per core was 262,144 (=643), and the total
problem size was 4,193,304. Table 1 shows the results of theperformance analyzer for Fujitsu FX10 [5]. Performance ofmultigrid process was evaluated. Performance of sliced ELL
is 23% larger than that of original ELL due to better
utilization of cache and effect of prefetching. Current versionof the code is not efficient, because the code stores two seriesofarrays
IAUnew3/AUnew3(3,N)
and
IAUnew6/AUnew6(6,N), where extra memory space isrequired.
Boundary Meshes
AUnew6(6,N)

introduced in recent works of parallel multigrid method on
peta-scale supercomputers with more than 105 cores, such as
IBM BlueGene/P, and Cray-XE6 systems [7,8]. Both of [7]and [8] adopted only flat MPI, and they dynamically reduce
the number of MPI processes in order to ensure that overheadof parallel communication does not dominate the computation.
Both of them define the minimum number of elements foreach process below which the communication costs start todominate. This critical value is defined as 2,000 in [7] and1,000 in [8]. In the present work, we define the switchinglevel according to the optimum one for CGA. If the optimum
switching level for CGA is levCGAopt, (=7 in the case of Fig.5),
the switching level to reduced number of processes (levhCGA) is
set to be levhCGA = levCGAopt -1, or =levCGAopt -2. If we consider
the case in Fig.5, levhCGA is set to 5 or 6. Moreover, we have to
specify the switching level to the final coarse grid solver(=levhCGA-CGA). In the present work, optimum value of levhCGA
and levhCGA-CGA are defined through empirical studies. Finally,
number of MPI processes after repartitioning (PErep) is
defined as 1/64 or 1/512 of the original number of MPI
processes. PErep is also decided empirically.
Fine
Level=1

Pure Internal
MeshesAUnew3(3,N)

Level=2

Level=mͲ3

Fig.11. Idea of sliced ELL format [12], 2 slices in the present work
TABLE 1: PERFORMANCE OF MULTIGRID PROCEDURE, MGCG SOLVER FOR
262,144 MESHES (=643) WITH A SINGLE NODE (16 CORES) OF FUJITSU FX10 BY
A PERFORMANCE ANALYZER [5], TOTAL PROBLEM SIZE IS 4,194,304

1.53u109

L2
Cache
Miss
1.67u1078

1.27u10

7

93.88

6.99

4.91u108

9.14u106

93.88

8.56

Instruction
CRS
Original
ELL (Fig.9)
Sliced ELL(Fig.10)

4.91u10SIMD
Operation
Ratio (%)
30.14

GFLOPS

Level=mͲ2

Coarsegridsolverona
singleMPIprocess(multiͲthreaded,furthermultigrid)

Coarse6.05

Fig.12. Procedures of hierarchical CGA (hCGA), where number of MPI
processes is reduced before the final coarse grid solver of CGA on a single
MPI process, In this case, number of MPI processes is reduced from 12 to 3 atlevel=m-3, and switched to the final coarse grid solver at level=m-2

G. hCGA for Parallel Communication
CGA (Coarse Grid Aggregation) [3] provides efficiency
and robustness of parallel multigrid method, as shown in Fig.5.
Problem of this approach is coarse grid problem is solved on asingle MPI process. Therefore, computational amount forcoarse grid solver increases as number of total cores increases,
as shown in Fig.7. This is more significant, if number of
threads for each MPI process is smaller, such as flat MPI, HB
4u4. In the present work, hierarchical version of CGA(hCGA) is introduced, as shown in Fig.12. In hCGA, numberof MPI processes is reduced and processes are repartitioned ina intermediate level before the final coarse grid solver on asingle MPI process. For example, number of MPI processes isreduced from 12 to 3 at level=m-3, and switched to the final
coarse grid solver at level=m-2, as shown in Fig.12. hCGA is
also expected to reduce communication overhead at coarser
levels of multigrid process. This type of approach is already

スーパーコンピューティングニュースLevel=mͲ3

IV. RESULTS
A. Overview
Optimization procedures for serial and parallel
communications in parallel multigrid described in theprevious chapter have been implemented for MGCG solver ofpGW3D-FVM [2,3,4]. Table 2 summarizes details of eachcase. Four types of parallel programming models (flat MPI,
HB 4u4, HB 8u2, and HB 16u1) are evaluated. Results of C3
and C4 are new in the present work.
TABLE 2: SUMMARY OF DETAILS OF EACH CASE

C0 [3]
C1 [4]
C2 [4]
C3
C4- 19 -

Format of Matrix Storage
CRS
ELL (Fig.9)
Sliced ELL (Fig.10)

CGA/hCGANO CGA (Fig.2)
CGA (Fig.3)
hCGA (Fig.12)

Vol.17,No.62015スーパーコンピューティングニュースTABLE 3: COMPARISON OF C3 (CGA) AND C4(HCGA) AT 4,096 NODES, FOR
FLAT MPI AND HB 4u4

Flat MPI
HB 4u4C3

levCGAopt
levhCGAopt
7

C4

6

C3

8

C4

6

levhCGA-CGAopt(PErep)3 (128 proc’s
8 nodes)
2 (256 proc’s,
64 nodes)Iter’s

sec.

64

13.2

61

8.22

59

8.08

56

7.97

15.0

Flat MPI:C3
HB 4x4:C3
HB 8x2:C3
HB 16x1:C3

sec.

12.5

10.0

7.5

5.0

100

1000

10000

100000

CORE#

Fig. 13. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), weak scaling: 262,144 (=643) meshes/core, max. total
problem size: 17,179,869,184 meshes, elapsed computation time for MGCG,
best results of sliced ELL with CGA (C3), switching to coarse grid solver atlevel=7 for CGA
20.00

HB 8x2:C0
HB 8x2:C1
HB 8x2:C2
HB 8x2:C3

sec.

15.00

10.00

5.00

0.00

100

1000

10000

100000

CORE#

Fig. 14. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), weak scaling: 262,144 (=643) meshes/core, max. total
problem size: 17,179,869,184 meshes, elapsed computation time for MGCG,
HB 8u2, C3 (sliced ELL) is 28% faster than C2 (ELL) at 4,096 nodes
15.0

10.0

Rest
Coarse Grid Solver
MPI_Allgather
MPI_Isend/Irecv/Allreducesec.

B. Weak Scaling
The performance of weak scaling was evaluated by using 8
to 4,096 nodes (65,536 cores) of the Fujitsu FX10. The
number of finite-volume meshes per core was 262,144 (=643);
therefore, the maximum total problem size was17,179,869,184 meshes.
Figure 13 compares four types of parallel programming
models for C3. Switching level to the coarse grid solver in
CGA is set to 7 for all of HB 4u4, HB 8u2 and HB 16u1,
which provides the best performance at 4,096 nodes for allcases. Each of HB 4u4, HB 8u2 and HB 16u1 provides
excellent scalability up to 4,096 nodes, although HB 8u2 isslightly faster than others. Flat MPI is much slower than otherhybrid parallel programming models with more than 104 coresbecause of overhead of coarse grid solver. Figure 14 provides
the performance of C0-C3 for HB 8u2. Switching level to thecoarse grid solver in CGA (C2 and C3) is set to 7.
Improvement of performance by sliced ELL over ELL(improvement from C2 to C3) at 4,096 nodes is 25%-31%,
and 28% for HB 8u2. Effect of sliced ELL is significant.
Figure 15 compares best results of C3 and C4 at 4,096
nodes. Effect of hCGA is not so significant except for flat
MPI where C4 is 1.61 times faster than C3. C4 is slightly
faster for HB 4u4, while C4 is slower than C3 for HB 8u2 andHB 16u1 due to the overhead of repartitioning for hCGA. 2digit number in each label after a colon on the transverse axisof Fig.15 indicates the number of iterations until convergence.
This means that hCGA can make changes on convergence.
Convergence of flat MPI and HB 4u4 is improved by hCGA.
Table 3 summarizes comparison between C3 (CGA) andC4 (hCGA) for flat MPI and HB 4u4. Optimum number of
MPI processes after repartitioning (PErep) for each case is 128
processes (8 nodes) for flat MPI, and 256 processes (64nodes) for HB 4u4, respectively. Optimum level for switching
to reduced number of MPI processes for hCGA (levhCGAopt) is
6 for both of flat MPI and HB 4u4, while levCGAopt for flat MPI
is 7, and that of HB 4u4 is 8. Critical mesh number per MPI
process for repartitioning is 1,000-2,000 in the existing works[7,8], but this value for the optimum case of flat MPI is 8 in
the present work, where number of meshes in each MPI
process is equal to 8 (=23) at lev=6 as shown in Table 3.
Overhead for coarse grid solver was significant for flat MPI
with more than 104 cores, as shown in Fig.7. At 4,096 nodes,
intermediate level of this coarse grid problem was solved by 8
nodes (128 MPI processes (=PErep)) for C4, instead of by asingle core for C2 and C3.
Figure 16 compares best results of C3 and C4 for fourparallel programming models. All of them except flat MPI(C3) are scalable up to 4,096 nodes, although HB 8u2 (C3) isslightly faster than others. Finally, Fig.17 provides theperformance of C3 and C4 for the four parallel programming
models at 64, 512, and 4,096 nodes. Generally speaking, C4
with hCGA shows significant improvement of performance offlat MPI, as number of nodes increases.

5.0

0.0

Flat MPI: Flat MPI: HB 4x4: HB 4x4: HB 8x2: HB 8x2: HB 16x1: HB 16x1:
C3:64
C4:61
C3:59
C4:56
C3:55
C4:56
C3:55
C4:56

Fig. 15. Performance of MGCG solver on Fujitsu FX10 using 4,096 nodes(65,536 cores), total problem size: 17,179,869,184 meshes, comparison of C3
and C4, effect of hCGA

- 20 -

Vol.17,No.62015boundaries. Generally speaking, hCGA improve the parallel
communications of MGCG solver. Effect of improvement is
15% for flat MPI, and less than 1% for HB 16u1. Finally,Fig.20 compares parallel communication overhead for strong
scaling computing at 4,096 nodes. Effect of improvement by
hCGA is more than that for weak scaling cases, where effectof improvement is 35% for flat MPI, and 8% for HB 16u1.

15.0

Flat MPI:C3
Flat MPI:C4HB 4x4:C4
HB 8x2:C3
10.0

HB 16x1:C3

120

7.5

5.0

100

1000

10000

Parallel Performance (%)

sec.

12.5

100000

CORE#

Fig. 16. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), weak scaling: 262,144 (=643) meshes/core, max. total
problem size: 17,179,869,184 meshes, elapsed computation time for MGCG,
best results of C3 and C4, and flat MPI (C3)15.0

C3, 64 nodes
C4, 64 nodes
C3, 512 nodes
C4, 512 nodes
C3, 4,096 nodes
C4, 4,096 nodes

80

Flat MPI:C3
HB 4x4:C3
HB 8x2:C3
HB 16x1:C3

Flat MPI:C4HB 4x4:C4
HB 8x2:C4
HB 16x1:C4

8192

65536

60
40
20
0

1024

CORE#

Fig. 18. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), strong scaling: 268,435,456 meshes (=1024u512u512),
parallel performance based on the performance of flat MPI (C3) with 8 nodes(128 cores)

sec.

10.0

100

5.0

5.00E-02

CGA-best (C3)hCGA-best (C4)4.00E-02
Flat MPI

HB 4x4

HB 8x2

HB 16x1

sec./iteration

0.0

Fig. 17. Performance of MGCG solver on Fujitsu FX10 using up to 4,096
nodes (65,536 cores), weak scaling: 262,144 (=643) meshes/core, max. total
problem size: 17,179,869,184 meshes, elapsed computation time for MGCG

2.00E-02
1.00E-02
0.00E+00

Flat MPI

HB 4x4

HB 8x2

HB 16x1

Fig. 19. Effect of hCGA on overhead of parallel communications in MGCG
solver, overhead for parallel communications per a single CG iteration, total
problem size: 17,179,869,184 meshes at 4,096 nodes
1.00E-02

CGA-best (C3)hCGA-best (C4)8.00E-03

sec./iteration

C. Strong Scaling
Finally, the performance of strong scaling was evaluatedfor a fixed size of problem with 268,435,456 meshes(=1024u512u512) using 8 to 4,096 nodes of Fujitsu FX10. At
4,096 nodes, the problem size per each core is only 4,096
meshes (=163). Figure 18 shows the performance of the
MGCG solver, and compares C3 (sliced ELL with CGA) andthe best case of C4 (sliced ELL with hCGA) for each parallel
programming model. The performance of flat MPI (C3) with
8 nodes (128 cores) is 100%, and the parallel performance is
50%–67% up to 512 nodes (8,192 cores). At 4,096 nodes, flat
MPI (C4) provides the best performance (14.6%).
Improvement of performance of flat MPI by hCGA (from C3
to C4) is significant, and the ratio of performance is 6.27
times at 4096 nodes. Improvement of performance by sliced
ELL over original ELL (from C2 to C3) at 4,096 nodes is
21.8% (flat MPI), 14.2% (HB 4u4), 15.3% (HB 8u2), and9.1% (HB 16u1), respectively.

3.00E-02

6.00E-03
4.00E-03
2.00E-03
0.00E+00

Flat MPI

HB 4x4

HB 8x2

HB 16x1

Fig. 20. Effect of hCGA on overhead of parallel communications in MGCG
solver, overhead for parallel communications per a single CG iteration, total
problem size: 268,435,456 meshes at 4,096 nodes

D. Improvement of Parallel Communication by hCGA
In this section, effect of hCGA on improvement of parallel
communication is evaluated for both of weak and strong cases.
V. SUMMARY AND FUTURE WORK
Figure 19 compares overhead for parallel communications per
The parallel multigrid and MGCG include both of seriala single CG iteration for weak scaling computation at 4,096 and parallel communication processes which are generally
nodes. Overhead for parallel communication includes, expensive. In the present work, new format for sparse matrix
MPI_Allreduce for dot products, MPI_Isend/Irecv/Waitall for storage based on sliced ELL is proposed for optimization ofcommunications at boundaries, and memory copies at

スーパーコンピューティングニュース- 21 -

Vol.17,No.62015serial communication on memories, and hierarchical coarse
grid aggregation (hCGA) is introduced for optimization ofparallel communication by message passing. The proposed
methods are implemented for pGW3D-FVM, and therobustness and performance of the code was evaluated using
up to 4,096 nodes (65,536 cores) of the Fujistu FX10 system.
The parallel MGCG solver using the sliced ELL format
provided performance improvement in both weak scaling(25%–31%) and strong scaling (9%–22%) compared to thecode using the original ELL format [4]. Moreover, hCGA
provided excellent performance improvement in both weakscaling (1.61 times) and strong scaling (6.27 times) for flat
MPI parallel programming model. hCGA was also effective
for improvement of parallel communications. Effect of sliced
ELL on serial communication was significant, while that ofhCGA on parallel communication was not so impressive
except for flat MPI cases. Because hCGA proved to be very
effective for reducing overhead of coarse grid solver, that will
also provide more significant effect on hybrid parallel
programming models with larger number of nodes. As was
mentioned in Chapter II, computational amount of coarse grid
solver for each core of flat MPI is 256 (=16u16) times as large
as that of HB 16u1. Therefore, hCGA is expected to be reallyeffective for HB 16u1 with more than 2.50u105 nodes(4.00u106 cores) of Fujitsu FX10, where the peak performanceis more than 60 PFLOPS.
In the present work, we focused on optimization ofcommunications with neighboring MPI processes at process
boundaries by MPI_Isend/Irecv/Waitall as parallel
communication, but overhead by collective communication(e.g. MPI_Allreduce) is more significant with more than
2.50u105 nodes. Figure 21 shows overhead by a single
MPI_Allreduce call for sending a single scalar variable for dot
product in the weak scaling case of MGCG solver, with up to
4,096 nodes of Fujitsu FX10. In order to reduce overhead bythis type collective communication, communication
avoiding/hiding type of algorithm combined with MPI
functions for non-blocking collective communications
supported in MPI 3.0 specification [13] is also effective.
7.00E-03

Flat MPI
HB 4x4
HB 8x2
HB 16x1

sec./MPI_Allreduce

6.00E-03
5.00E-03
4.00E-03
3.00E-03
2.00E-03
1.00E-03
0.00E+00

100

1000

10000

100000MPI Process #Fig. 21. Overhead by MPI_Allreduce for dot product (a single scalar variable)of MGCG solver on Fujitsu FX10 using up to 4,096 nodes (65,536 cores),weak scaling: 262,144 (=643) meshes/core, max. total problem size:
17,179,869,184 meshes, best results of C3 and C4More efficient storage of sliced ELL arrays andoptimization of hCGA function are the critical issues for

スーパーコンピューティングニュースfuture work. Especially, improvement of communication
procedures for repartitioning of MPI processes is important.
CGA and hCGA include a various types of parameters, andthe optimum values of those were derived through empiricalstudies in the present work. Development of methods forautomatic selection of these parameters [14] is also aninteresting technical issue for future work. Optimum
parameters can be estimated based on calculation ofcomputational amounts, performance models, parameters of
hardware, and some measured performance of the system. But
it is not so straightforward. Because some of these parameters
also make effects on convergence, construction of suchmethods for automatic selection is really challenging.
ACKNOWLEDGMENT
This work is supported by Core Research for Evolutional
Science and Technology (CREST), the Japan Science andTechnology Agency (JST), Japan. The computational
resource of Fujitsu FX10 was awarded by the “Large-scale
HPC Challenge” Project, Information Technology Center, theUniversity of Tokyo.
REFERENCES[1] HPCG: High Performance Conjugate Gradients:
https://software.sandia.gov/hpcg/[2] Nakajima, K., New strategy for coarse grid solvers in parallel multigrid
methods using OpenMP/MPI hybrid programming models”, ACM
Proceedings of the 2012 International Workshop on Programming
Models & Applications for Multi/Manycores (2012)[3] Nakajima, K., OpenMP/MPI Hybrid Parallel Multigrid Method on
Fujitsu FX10 Supercomputer System, IEEE Proceedings of 2012
International Conference on Cluster Computing Workshops, IEEE
Digital Library: 10.1109/ClusterW.2012.35 (2012) 199-206[4] Nakajima, K., Large-scale Simulations of 3D Groundwater Flow using
Parallel Geometric Multigrid Method, Procedia Computer Science 18(2013) 1265-1274[5] Information Technology Center, The University of Tokyo:
http://www.cc.u-tokyo.ac.jp/[6] Baker, A., T. Gamblin, M. Schultz, U. Yang, Challenge of Scaling
Algebraic Multigrid across Modern Multicore Architectures,
Proceedings of the 2011 IEEE International Parallel & Distributed
Processing Symposium (IPDPS’11) (2011) 275-286[7] Lin, P.T., Improving multigrid performance for unstructured meshdrift-diffusion simulations on 147,000 cores, International Journal for
Numerical Methods in Engineering 91 (2012) 971-989[8] Sundar, H., G. Biros, C. Burstedde, J. Rudi, O. Ghattas, and G. Stadler,
Parallel Geometric-Algebraic Multigrid on Unstructured Forests ofOctrees, ACM/IEEE Proceedings of the 2012 International Conference
for High Performance Computing, Networking, Storage and Analysis(SC12) (2012)[9] Fujitsu: http://www.fujitsu.com/[10] Deutsch, C.V. and A.G. Journel, GSLIB Geostatistical Software
Library and User’s Guide, Second Edition. Oxford University Press(1998)[11] Smith, B., P. BjIrstad, and W. Gropp, Domain Decomposition,
Parallel Multilevel Methods for Elliptic Partial Differential Equations,Cambridge Press (1996)[12] Monakov, A., A. Lokhmotov, and A. Avetisyan, Automatically tuning
sparse matrix-vector multiplication for GPU architectures, Lecture
Notes in Computer Science 5952 (2010) 112-125[13] Ghysels, P. and W. Vanroose, Hiding global synchronization latency inthe preconditioned Conjugate Gradient algorithm, Parallel Computing
40-7 (2014) 224-238[14] Nakajima, K., Automatic Tuning of Parallel Multigrid Solvers using
OpenMP/MPI Hybrid Parallel Programming Models, Lecture Notes in
Computer Science 7851 (2013) 435-450- 22 -

Vol.17,No.62015