1jh210026-NAH

Innovative Multigrid Methods II
Akihiro Fujii（Kogakuin University）
Abstract
Multigrid method is known as one of the most eﬃcient and scalable methods using thesupercomputers. It is applicable to partial diﬀerential equations with ﬁne and coarse grids,
linear matrix problems, and parallel time integration. Our project studies multigrid related
themes in 3-year plan. This report is for the second year. Although we continue studying
the area, we’ve conducted researches in 2021, including the topics of precision optimization with SELL-C-σ sparse matrix format, computation and communicaiton overlapping insparse matrix vector product, weak scaling evaluation of AMG solver, new parareal method
for explicit time-marching applications, and scalability enhancement method for Parallel inSpace Time solvers.

1 Basic Information
1.1 Collaborating JHPCN Centers
Tokyo, Nagoya
1.2 Research Areas- Very large-scale numerical computation
1.3 Roles of Project Members
Akihiro Fujii3 :AMG, PinST
Kengo Nakajima2,4 : (Co-PI) Application,
GMG, PinST, AMG
Matthias Bolten8 : (Co-PI) GMG, PinST,
AMG
Masatoshi Kawai2 : GMG, AMG
Akihiro Ida13 : GMG, AMG
Gerhard Wellein9 : GMG, AMG
Christie Alappat9 : GMG, AMG
Martin Schreiber11 : GMG, AMG
Tetsuya Hoshino2 : GMG, AMG
Satoshi Ohshima6 : GMG, AMG
Toshihiro Hanawa2 : GMG, AMG
Osni Marques10 : GMG, AMG
Kenji Ono5 : PinST
Mikio Iizuka5 : PinST
Takeshi Iwashita1 : AMG, PinST

Yasuhito Takahashi7 : PinST
Robert Speck12 : PinST
Atsuhiro Miyagi14 : PinST
Teruo Tanaka3 :AMG, PinST
Alexander T.Magro2 : GMG
Ryo Yoda2 : PinST, AMG
Yen-Chen Chen2 : PinSTGayatri Caklovic12 :PinST
1: Hokkaido U., 2: U. Tokyo, 3: Kogakuin U.,
4: RIKEN R-CCS, 5: Kyushu U., 6: Nagoya U., 7:Doshisha U., 8: U. Wuppertal*, 9: FAU*, 10: LBNL,USA, 11: Technical U. of Munich*, 12: Juelich Supercomputing Centre*, 13: JAMSTEC, 14: Taisei
Corp. *:Germany2 Purpose and signiﬁcance ofResearchMultigrid method is a promising approach forlarge-scale computing in exa-scale era. We
develop robust and eﬃcient parallel multigrid methods for both of geometric multigrid (GMG) and algebraic multigrid(AMG),
focusing on robust and eﬃcient smoothers,
and hierarchical methods which are proposedand developed by ourselves. We have original PinST methods. Developed methods will

Final Report for JHPCN Joint Research of FY 2021

2

be implemented as a numerical library andit will be applied to various types of applications. It is expected to provide outstanding
performance for large-scale real-world applications.

1012 unknowns in [AMG 1]. As for PinST,
we studied parareal method for explicit timemarching scheme in [PinST exp 1]. Performance enhancement technique adjusting thebalance of parallelism in space and time directions was proposed in [PinST scl 1].

3 Signiﬁcance as JHPCN Joint
Research Project
Multigrid method is scalable and used inmany ﬁelds. It is known as one of the mosteﬃcient linear solvers. It can also be applied to parallel time integration problems,
which exploits parallelism in time dimension. Our research project has original codes
and algorithms as written in Research purpose. Therefore, research papers and codesfrom the project will enhance the eﬃciencyof the multigrid solver, and will help many
researchers exploit parallelism in time direction.
Our research focuses on hierarchical algorithms and their performance on supercomputers. Thus, availability of supercomputers with diﬀerent kinds of architectures
helps us verify the codes we are developing. In addition, a JHPCN joint research
project oﬀers collaborative research opportunity with JHPCN members who have expertise knowledge in various application ﬁelds.
Our project members include international
experts in Germany, US, and Japan on multigrid methods and PinST. We are sure that
this JHPCN joint research project promotes
the international collaborative activity with
JHPCN members.

4 Outline of Research Achievements
up to FY2021
This project is set as a 3-year project, andthis report is for the second year. Main presentation and publication of the second year
is listed in section 7. For GMG, we have conducted precision optimization with SELL-Cσ in [GMG 1], and computation and communication overlapping in [GMG 4]. For AMG,
we have checked our AMG solver’s scaling
performance with large sized problems with

5 Details of FY2021 ResearchAchievements
In the 2nd year of our continuing project, weare doing researches including the following
items.
1. GMG: Precision optimization withSELL-C-σ
2. AMG: Weak scaling performance test on
Supercomputer Flow and Fugaku
3. PinST: New explicit time-marching application
4. PinST: Scalability enhancement forMGRIT
5. Comm: Fast halo communication directly using RDMA interface
Subsections introduce the researches.
5.1 GMG: Parallel multigrid using SELL-C-σwith FP64/32, and communication computation overlapping technique
A paper on parallel multigrid method using SELL-C-σ with FP64/FP32 has beenpresented at iWAPT workshop in IPDPS
2021 [GMG 1], and there were many usefulcomments from audience for future works,
such as eﬀects of precision on MPI communications. The developed method in [GMG 1]
has been applied to Wisteria/BDEC-01(Odyssey) with A64FX, and performancehas been evaluated using up to 4,096
nodes [GMG 2,3]. Conﬁguration of network
topology aﬀects the performance of all-to-alltype communications, if the number of nodes
is more than 2,048.
In [GMG 4], we studied communicationcomputation overlapping. Communication
overhead is a critical issue when executing sparse matrix solvers on large-scale massively parallel supercomputers. In the previous work, we introduced communicationcomputation overlapping with dynamic loop

25.0

30.0
OFP-L(a)Final Report for JHPCN
Joint Research of FY 2021
OFP-M

20.0

OFP-S

Improvement (%)

Improvement (%)

25.0

(a)

20.0

15.0

10.0

Static Dyn: 200 300 400 500 600 700 800 900 1000
100

30.0

15.0

10.0

5.0

25.0

(b)

ODY-L
ODY-M
ODY-S

20.0

Improvement (%)

Improvement (%)

3(5.0)

0.0

25.0

PJ-M
CRS-M
ELL-M

0.05.0(b)

PJ-L
CRS-L
ELL-L

20.0

15.0

10.0

15.0

PJ-L
PJ-M
CRS-L
CRS-M
ELL-L
ELL-M

10.0

5.0

0.0
5.0(5.0)
0.0
Static Dyn: 200 300 400 500 600 700 800 900 1000
100

Fig. 1 Improvement by Static and Dynamic (Dyn) CC-Overlapping over “original” implementations at 2,048 nodes (a)OFP with IHK/McKernel, (b) Odyssey

Fig. 2 Improvement by Static and Dynamic CC-Overlapping over “Original”
implementation at 1,024 nodes (a) OFP
with IHK/McKernel, (b) Odyssey

5.2
scheduling of OpenMP to the sparse matrixvector multiplication (SpMV) process of aparallel iterative solver by Conjugate Gradient (CG) method in a parallel ﬁnite element application (GeoFEM/Cube) on multicore and manycore clusters. In the presentwork, ﬁrst, we re-evaluated the method
on Wisteria/BDEC-01 (Odyssey) (Fujitsu
PRIMEHPC FX1000 with A64FX), and asigniﬁcant performance improvement of 2530% for parallel iterative solver at 2,048
nodes (98,304 cores) was obtained. Moreover, we proposed a new reordering method
for communication-computation overlappingin ICCG solvers for a parallel ﬁnite volume application (Poisson3D/Dist), and attained 5-12% improvement at 1,024 nodes ofOdyssey.

AMG: Weak scaling performance test on
Supercomputer Flow and Fugaku
We proposed light-weight process aggregation method using parallel C/F splitting algorithm (CLJP). From the numerical tests,
it turned out that our CLJP-based method
was better than Coarse grid aggregation withParMETIS library especially in the computing environment with more than 100000 processes. In addition to the new process aggregation method, we improved our implementation for large sized coarse level matrix
whose size is over the range of 32-bit integertype.
We evaluated the solver with this CLJPbased coarse grid aggregation on Supercomputer Fugaku using up to about 50000 nodes.
Here, we introduce the weak scaling evaluation with Supercomputer Fugaku. Numericaltest condition is listed in Table 1. Since more

Final Report for JHPCN Joint Research of FY 2021

4

Table 1 Weak scaling test setting

Domain size ofa process
Process/thread
allocation on 1 node
CRS matrix sizeof 1 node
Number of nodes
Max.number ofcores and processes
Max. domain size

288 × 144 × 144
4 MPI processes
12 Open-MP thrds
7.8GB27 ele. per row
from 768 to49152 nodes
2.35 × 106 cores
1.9 × 105 processes
1.17 × 1012 DOF

than 50000 nodes are used in the test, 4 processes are assigned to each node. The CPU
has 4 CMGs (Core Memory Groups) and 4
MPI-processes per 1 node setting is recommended in general. The problem matrix size
was set as about 8GB per node, which occupies a 1/4 each node 32GB memory.
The solver was AMG preconditioned
CG method with multi-color Gauss Seidelsmoother. Thus it was executed in MPIOpenMP hybrid model. The convergence criterion was the relative residual 2-norm lessthan 1.0 × 10−6 .
Fig. 3 and Fig. 4 show total execution
time and AMG-CG one iteration time, respectively. The horizontal axis of these ﬁgures indicates the problem DOFs in logarithmic scale. The vertical axis indicates thetime in seconds. Fig. 3 plots total execution
time and multi-level creation (Setup) time.
Numbers in the ﬁgure correspond to the iteration number for convergence. Total execution time requires only 1.31 times longer
for 64 times larger problem. Iteration number also increases from 30 to 37. As a whole,it demonstrates a good weak scaling performance. Fig. 4 plots AMG-CG one iteration
time. From smallest to the largest problem
sizes, the time diﬀerence was only 7 %.
In the largest problem, the number of unknowns at the second coarse level is 4.32 ×1010 which is over the maximum value of 32

Fig. 3 Weak scaling test on Fugaku

Fig. 4 AMG-CG 1 iteration timebit integer. For such problems, the solver
creates distributed matrices with 32-bit integer, and reaches convergence. Hence, it isassumed to be applicable to larger problems,
as long as the block row matrix allocated to1 process is represented in the range of 32-bit
integer type.
5.3 PinST: explicit time-marching application
Parallel-in-Space/Time (PinST) method forexplicit time-marching application has beendeveloped based on parareal method. The
method has been applied to 1D convectiondiﬀusion equations, and 2D compressible viscous ﬂow problems with Navier-Stokes equations. Proposed method provides better performance and stability compared to existingmethods, such as parareal, and MGRIT. This

Final Report for JHPCN Joint Research of FY 2021

Fig. 5 Execution time comparison of
spatial parallelization to parallel-inspace/time for one-dimensional advection
example.

is one of the ﬁrst example where PinST is
successful for practical applications with explicit time-marching. The proposed method
and results were presented at domestic andinternational meetings [PinST exp 2,3], anda paper submitted to ScalA21 workshop
in SC21 has been accepted for presentation [PinST exp 1]. While the method is applied to subsonic ﬂow currently, method forsupersonic cases using adaptive mesh reﬁnement is under development. Fig. 5 introduces
the performance comparison among Parallelin-space, Parallel-in-time, and Parallel-inspace/time for one-dimensional advection example.
5.4 PinST: MGRIT scaling enhancement
In the previous work, we apply coarse-grid
agglomeration (CGA) to MGRIT to enhance
scalability, which gradually shrinks the number of active processes on a coarse grid.
While this technique achieved well-scaling
performance, it makes more processors idle
than usual spatial multigrid methods. Therefore, if these idle processors can be utilized
somehow, it will lead to further performance
improvement.
These researches [PinST scl 1,2,4] consider
a method to assign these idle processors to
the spatial redistributed domain on coarse
levels. Then, temporal parallelism decreases
by CGA, and spatial parallelism increases
with each level. This method aims to accel-5

erate coarse-level spatial solvers.
Numerical experiments were conducted forthree-dimensional diﬀusion problems on astructured grid. We observed the well-scaling
performance of up to 512 nodes on OFP andthat the proposed method achieves about
1.4 times faster than pure MGRIT. We conﬁrmed the eﬀectiveness of the use of idle processors by CGA.
In addition, the research [PinST scl 3] investigates the optimized coarse grid operator
for time-dependent Stokes and Oseen problems.
5.5 Fast halo Communication
Supercomputers such as Fugaku and FX1000
oﬀer users RDMA communication interface.
We investigated fast halo communication using this interface in [Comm 1] and [Comm 2].
One-sided communications such as RDMA
write tends to need explicit synchronization in general. We eliminated this explicit
synchronization by preparing double buﬀers,
and evaluated the eﬃciency of the halo communication in typical communication setting.

6 Progress during FY2021 and Future
Prospects
Three areas of GMG, AMG and PinST are
set in our project. We have made progressin all areas, and the progress is going well.
In the next ﬁscal year, we will organize theresearch results to publish papers and codes.

7 List of Publications andPresentations
Journal Papers (Refereed)
Proceedings of International Conferences (Refereed)[PinST scl 1] Yoda, R., Bolten, M.(+),
Nakajima, K., Fujii, A., Assignment ofidle processors to spatial redistributed
domains on coarse levels in multigrid reduction in time. In HPC Asia 22: The
International Conference on High Performance Computing in Asia-Paciﬁc Region, Jan 12–14, 2022, Kobe, Japan &Online. ACM, New York, NY, USA (inFinal Report for JHPCN Joint Research of FY 20216press)[PinST exp 1] Chen, Y.-C., Nakajima, K.,
Optimized Cascadic Multigrid Parareal
Method for Explicit Time-Marching
Scheme, IEEE Proceedings of ScalA21:
12th Workshop on Latest Advances in
Scalable Algorithms for Large-Scale Systems, in conjunction with SC21, November 2021 (in press)[GMG 1] Nakajima, K.,
Ogita,T.,
Masatoshi, K., Eﬃcient Parallel Multigrid Methods on Manycore Clusters
with Double/Single Precision Computing, IEEE Proceedings of iWAPT 2021
in conjunction with IPDPS 2021, May2021[GMG 4] Horikoshi,
M.,Geroﬁ,
B.,
Ishikawa, Y., Nakajima, K., Exploring
Communication-Computation
Overlap in Parallel Iterative Solvers onManycore CPUs using Asynchronous
Progress Control, IXPUG Workshop inconjunction with HPC Asia 2022, 2022
[GMG 5] Kawai, M.,
Nakajima, K.,
Low/Adaptive Precision Computation in Preconditioned Iterative Solvers
for Ill-Conditioned Problems, TheInternational Conference on High Performance Computing in Asia-PaciﬁcRegion (HPC Asia 2022), 2022

Time Integration, August 2021
[PinST scl 2] Yoda, R., Bolten, M.(+),
Nakajima, K., Fujii, A., Spatial redistribution on the temporal coarse levelfor Multigrid Reduction in Time, PinT
2021: 10th Workshop on Parallel-inTime Integration, August 2021
[PinST scl 3] 依 田 凌, 中島研吾, 藤 井 昭 宏,
MGRIT の粗格子演算子に対する最適化手法とその線形時間発展 Stokes ・ Oseen
問題への適用, 日本応用数理学会「行列・Proceedings of International Conferences (Nonrefereed)[GMG 2] 中島 研吾，河合 直聡，
Wisteria/BDEC-01(Odyssey) に お ける並列多重格子法ソルバーの開発と性能評
価，日本応用数理学会年会 2021，先進的環境における数値計算と関連 HPC 技術，
2021 年 9 月[AMG 1] 藤井 昭宏，田中 輝雄，中島 研吾，
SA-AMG 法における軽量な粗格子集約手法と富岳上でのウィークスケーリング性能評価，情報処理学会研究報告（2021-HPC181）(第 181 回 HPC 研究会)（オンライン，2021 年 9 月 27 日）[PinST exp 2] Chen, Y.-C., Nakajima, K.,
An Eﬃcient Parallel-in-Time Method forExplicit Time-Marching Schemes, PinT
2021: 10th Workshop on Parallel-in-固有値問題の解法とその応用」研究部会第32 回研究会. (オンライン, 2021 年 12 月10 日 )[PinST scl 4] Yoda, R., Bolten, M.(+),
Nakajima, K., Fujii, A., Assignment ofIdle Processors to Spatial Redistributed
Domains for Multigrid Reduction inTime, SIAM Conference on Parallel Processing for Scientiﬁc Computing (PP22),
February 23–26, 2022, Virtual Conference, Originally scheduled in Seattle,
Washington, U.S.
[GMG 3] 中島 研吾，河合 直聡，
Wisteria/BDEC-01(Odyssey) に お ける大規模前処理付き反復法ソルバーの性能評価，日本応用数理学会 「行列・固有値問題の解法とその応用」研究部会（MEPA），
2021 年並列／ 分 散／ 協 調処理 に 関す る
サマー・ワークショップ (SWoPP2021)，
2021 年 7 月[PinST exp 3] Chen, Y.-C.,Nakajima, K.,
An Eﬀective Parallel-in-Time Method
for Explicit Time-Marching Schemes,IPSJ，IPSJ SIG Technical Report（2021-HPC-180-2），2021[PinST exp 4] Chen, Y.-C., Nakajima, K.,
Cascadic Parareal Method for Explicit
Time-Marching Schemes, HPC Asia
2022 (3rd Prize, HPC Asia 2022 Best
Student Poster Award in Memory of Hiroshi Nakashima), 2022
[Comm 1] Yoshimoto, K., Fujii, A., Tanaka,T., Fast Adjacent Communication withRDMA and Double Buﬀering, The 20th
International Symposium on Advanced
Technology (ISAT-20), Lightning Talk
session, Nov.23-24, 2021.
[Comm 2] Yoshimoto, K., Fujii, A., Tanaka,Final Report for JHPCN Joint Research of FY 2021

T., RDMA with Double Buﬀering for
Adjacent Communication, The International Conference on High Performance
Computing in Asia-Paciﬁc Region (HPCAsia 2022), poster, Jan. 2022.
Published library and relating data
Other (patents, press releases, books and soon)

7