---
title: 第5回 配信講義　計算科学技術特論A （2023）
tags: 
author: [R-CCS　計算科学研究推進室](https://convert.docswell.com/user/2300203199)
site: [Docswell](https://www.docswell.com/)
thumbnail: https://bcdn.docswell.com/page/972GM6VZ7R.jpg?width=480
description: 第5回  5月18日　プログラム高速化実例と大規模学習への展開 プログラム高速化の基礎知識、並列化プログラミング（MPI、OpenMP）の基礎知識、およびプログラム高速化の応用事例の座学を通して、 計算科学で必要な高性能計算技術の基礎の習得を目指す。
published: May 15, 23
canonical: https://convert.docswell.com/s/2300203199/ZP9MYE-2023-05-15-155045
---
# Page. 1

![Page Image](https://bcdn.docswell.com/page/972GM6VZ7R.jpg)

内容に関する質問は
katagiri@cc.nagoya-u.ac.jp
まで
第５回
プログラム高速化実例と
大規模学習への展開
名古屋大学情報基盤センター 片桐孝洋
1
2023年度
計算科学技術特論A


# Page. 2

![Page Image](https://bcdn.docswell.com/page/DJYVX93NJM.jpg)

講義日程と内容について

２０２３年度 計算科学技術特論Ａ（木曜：１３：００－１４：３０ ）





2
第１回：プログラム高速化の基 礎、2023年4月13日
 イントロダクション、ループアンローリング、キャッシュブロック化、
数値計算ライブラリの利用、その他
第２回：MPIの基礎、2023年4月20日
 並列処理の基礎、MPIインターフェース、MPI通信の種類、その他
第３回：OpenMPの基礎、2023年4月27日
 OpenMPの基礎、利用方法、その他
第４回：Hybrid並列化技法(MPIとOpenMPの応用)、2023年5月11日
 背景、Hybrid並列化の適用事例、利用上の注意、その他
第５回：プログラム高速化実例と大規模学習への展開、2023年5月18日
 プログラムの性能ボトルネック に関する考えかた（I/O、単体性能
(演算機ネック、メモリネック)、並列性能(バランス))、性能プロファイル、
機械学習におけるHPC、ほか
2023年度 計算科学技術特論A


# Page. 3

![Page Image](https://bcdn.docswell.com/page/V7NM1LVVE8.jpg)

性能チューニングの応用
3
2023年度
計算科学技術特論A


# Page. 4

![Page Image](https://bcdn.docswell.com/page/YJ9N546VJ3.jpg)

性能チューニングに関する総論（その１）
 コンパイラを過信しない
 書き方が悪いと、自動並列化だけでなく、
逐次最適化もできない！
 １ループ中で書いてある＜式＞がとても多い時
 レジスタが足りなくなって、メモリにデータを
吐き出すコードを生成するので、性能低下する
加えて、式が複雑で、コンパイラの解析ができない
要因が増える
⇒後述の、手による「ループ分割」が必要になる
4
2023年度
計算科学技術特論A


# Page. 5

![Page Image](https://bcdn.docswell.com/page/GJ8KYQZ4ED.jpg)

性能チューニングに関する総論（その２）

コンパイラを過信しない（つづき）
 自動並列化は＜特に＞過信しない
 ループ並列性がない逐次コードは並列化できない
 書き方が悪いと、原理的に並列化できるループも、
自動並列化できない
ループの構造（開始値、終了値が明確か、など）
言語的特徴から生じる問題もある
 C言語では、並列化したいループがある関数コール時の
引数にデータ依存があると判断されると、並列化できない。
自動並列化コンパイラ
！＝逐次コードを自動で並列コードに書き直すコンパイラ
5
2023年度
計算科学技術特論A


# Page. 6

![Page Image](https://bcdn.docswell.com/page/LJLZ4XDV7R.jpg)

性能チューニングに関する総論（その２）

コンパイラを過信しない（つづき）

例）foo (A, B, C); ←一般にA、B、Cは同一配列
で引渡される可能性があるため、A、B、C間は
依存があるとデフォルトで仮定
⇒ディレクティブ、コンパイラオプション指定で対応
int foo(double A[N][N], double B[N][N], double C[N][N]) {
int i, j, k;
for (i=0; i&lt;N; i++) {
for (j=0; j&lt;N; j++) {
for (k=0; k&lt;N; k++) {
C[i][j] += A[i][k] * B[k][j]; } } }
}
6
2023年度
計算科学技術特論A


# Page. 7

![Page Image](https://bcdn.docswell.com/page/47MVKLPN7W.jpg)

性能チューニングに関する総論（その３）

コンパイラを過信しない（つづき）
 スレッド数の増加

低スレッド並列（２～４スレッド）向きのコードと、高スレッド並列
（１６スレッドを超える）向きコードは、まったく異なる
コンパイラは、実行前にユーザが使うスレッド数を
知ることが出来ない
平均的なスレッド数を仮定、まあまあな
性能のコードを生成する
 並列数が増加すると、ループ長が短くなることで、
ループ並列性が無くなる
⇒後述の、手による「ループ融合」が必要になる

7
2023年度
計算科学技術特論A


# Page. 8

![Page Image](https://bcdn.docswell.com/page/P7R2XK14J9.jpg)

性能チューニングに関する総論（その４）

コンパイラを過信しない（つづき）



8
あるベンダ提供のコンパイラで最適化できたとしても、
別のベンダ提供のコンパイラで最適化できる保証はない
 例）富士通Fortranコンパイラ と インテルFortranコンパイラ
同一ベンダのコンパイラでも新規ハードでは、
同一コードでも同じ品質で最適化できる保証がない
従来からあるコード（レガシーコード）で、ハードウェア、および、
ソフトウェア環境が変わっても、高い性能を保つこと（性能可搬性）は、
HPC分野でホットな研究テーマ
 「ソフトウェア自動チューニング」の研究分野
 ソフトウェア性能工学 (Software Performance Engineering, SPE)
 ソフトウェア開発コストを低く保つ、チューニングの枠組み
 コード自動生成技術
 性能モデリング、パラメタ最適化、機械学習、等の技術が必要
2023年度
計算科学技術特論A


# Page. 9

![Page Image](https://bcdn.docswell.com/page/PJXD6LYZEX.jpg)

性能チューニングに関する総論（その５）

自分のコードのホットスポット（重い部分）を
認識せよ
 自分のコードのうち、どの部分が重いのか、
実測により確認せよ
1. 演算時間ボトルネック（演算時間が多い）
2. 通信時間ボトルネック（通信時間が多い）
3. I/Oボトルネック（I/O時間が多い）
9
2023年度
計算科学技術特論A


# Page. 10

![Page Image](https://bcdn.docswell.com/page/3JKX1L3VED.jpg)

性能チューニングに関する総論（その６）

自分のコードのホットスポット（重い部分）を認識せよ
 計算量など、机上評価はあてにならない
 実性能は計算機環境や実行条件に依存
思わぬところに ホットスポット（重い部分） あり
チューニング状況に応じホットスポットは変わる
 計算量が多くても、問題サイズが小さく、キャッシュに
のる場合は、演算時間が占める割合は少ない
（かもしれない）
 通信量が少なくても、通信＜回数＞が多いと、
通信レイテンシ律速
 I/O量が少なくても、I/Oハードウェアが貧弱、
実行時に偶発的にI/O性能が劣化すると、I/O律速
10
2023年度
計算科学技術特論A


# Page. 11

![Page Image](https://bcdn.docswell.com/page/LE35Q3LQJ5.jpg)

状況に応じて変化していくホットスポット

最初は演算律速
演算時間

通信時間
I/O
時
間
演算チューニングをすると、次は通信律速に
演算
時間
11
通信時間
I/O
時
間
2023年度
計算科学技術特論A


# Page. 12

![Page Image](https://bcdn.docswell.com/page/8EDD94NWEG.jpg)

ホットスポット判明後の最適化方針の一例

演算ボトルネックの場合（順番は検討する優先度）
1. コンパイラオプションの変更


プリフェッチ、ソフトウェア・パイプライニング強化オプション、など
アンローリング、タイリング（ブロック化）のディレクティブ追加、など
アルゴリズムを変更し、計算量が少ないものを採用
アルゴリズムを変更し、キャッシュ最適化向きのもの
を採用
2.
3.

「ブロック化アルゴリズム」の採用
コンパイラが自動で行わないコードチューニングを
手で行う
4.


12
アンローリング、ループ分割、ループ融合、など
高速化（連続アクセス）に向くデータ構造を採用
2023年度
計算科学技術特論A


# Page. 13

![Page Image](https://bcdn.docswell.com/page/V7PXZM5XE8.jpg)

ホットスポット判明後の最適化方針の一例

通信ボトルネックの場合
 通信レイテンシが主要因（通信回数が多い）
1. こま切れの通信をまとめて送る
（通信のベクトル化）
2. 冗長計算による通信回数の削減
3. 非同期通信による通信の隠ぺい
 通信量が主要因（１回当たり通信データが多い）
1. 冗長計算による通信量の削減
2. より高速な通信実装方式の採用
（Remote Direct Memory Access (RDMA) など）
3. 非同期通信による通信の隠ぺい
13
2023年度
計算科学技術特論A


# Page. 14

![Page Image](https://bcdn.docswell.com/page/2JVLW9P37Q.jpg)

ホットスポット判明後の最適化方針の一例

I/Oボトルネックの場合
1.


2.
3.

4.


14
「富岳」では、１ノード48MPIプロセスが走る。
ピュアMPIで、１MPIあたり１つファイルアクセス
をすると、IO性能が極端に劣化する
→１ノード当たり30MPIプロセス以下にする
高速なファイルシステムを使う
ファイルステージングの利用
高速ファイルシステム（ＳＳＤ、バーストバッファなど）の利用
データを間引き、I/O量を削減する
OSシステムパラメタの変更
I/Oストライプサイズの変更
 大規模データサイズを１回I/Oする場合は、ストライプサイズを大きくする
より高速なI/O方式を採用する
ファイル書き出しは、MPIプロセスごとに別名を付け、同時にI/O出力する
実装であることが多い
高速なファイルI/O（Parallel I/O、MPI-IOなど）を使う
 複数のファイルを１つに見せることができる
 ただしMPI-IO等が高性能でないと意味が無い
2023年度
計算科学技術特論A


# Page. 15

![Page Image](https://bcdn.docswell.com/page/5EGXDZPYJL.jpg)

ホットスポットをどのようにして知るのか
1.
2.
プログラム中にタイマを設定して調べる
性能プロファイラを利用する
 演算ボトルネック



通信ボトルネック


プロファイラの基本機能により調査可能
 例）富士通 基本プロファイラ、など
I/Oボトルネック


15
プロファイラの基本機能により調査可能
ループごとの詳細プロファイルにより、ハードウェア性能
（キャッシュヒット率など）を調査可能
 例）富士通 基本プロファイラ、など
一般にあまり提供されていない
スパコンベンダによっては専用プロファイラを提供している
 例）HPE（Cray）社のプロファイラ（CrayPat Performance Analysis Tool）
2023年度
計算科学技術特論A


# Page. 16

![Page Image](https://bcdn.docswell.com/page/4JQ81L267P.jpg)

その他の注意

I/Oを行うため、プロセス0にデータを集積し、
プロセス0のみがI/Oをするプログラム


データ集積のために、MPI_AllgatherV関数などが使われる
I/Oのための通信時間が占める割合が大きくなる
 ノード数が増えるほど、上記のI/O時間の割合は大きくなる
⇒超並列向きではない実装
I/Oは、プロセスごとに並列に行うほうが良い



16
ただし、プロセスごとに分散されて生成されるファイルの
扱いが問題になる
できるだけ、MPI-IOや、その他のシステムソフトウェア提供の
機能を使い、プロセスごとにファイルを見せない実装がよい
2023年度
計算科学技術特論A


# Page. 17

![Page Image](https://bcdn.docswell.com/page/K742RDLMJ1.jpg)

性能プロファイリング
17
2023年度
計算科学技術特論A


# Page. 18

![Page Image](https://bcdn.docswell.com/page/LJ1MQZLYEG.jpg)

性能プロファイリングの重要性

プログラムにおいて、どの箇所（手続き（関数））に時間
がかかっているか調べないと、チューニングを行っても
効果がない



手続きA:100秒、手続きB：10秒、手続きC：１秒、全体：１１１秒
手続きAは全体時間の90％なので、これをチューニングすべき
性能プロファイルを行うには、一般的には、スパコン提供
メーカが提供しているプロファイラを使うとよい

多くは、コンパイラと連携している
1.
2.
3.
4.
18
コンパイラオプションで指定し、実行可能コードを生成
実行可能コードを実行
性能プロファイルのためのファイル（ログファイル）が作成される
専用のコマンドを実行する
2023年度
計算科学技術特論A


# Page. 19

![Page Image](https://bcdn.docswell.com/page/GJW3M9V1J2.jpg)

性能プロファイラでわかること
性能プロファイラツールに大きく依存
 ノード内性能







全体実行時間に占める、各手続き（関数）の割合
MFLOPS（GFLOPS）値
キャッシュヒット率
スレッド並列化の効率（負荷バランス）
I/O時間が占める割合
ノード間性能

19
MPIなどの通信パターン、通信量、通信回数
（多くは専用のGUIで見る）
2023年度
計算科学技術特論A


# Page. 20

![Page Image](https://bcdn.docswell.com/page/4EZVQ94XJ3.jpg)

性能プロファイラ例





富士通コンパイラには、性能プロファイラ機能がある
富士通コンパイラでコンパイル後、実行コマンドで指定し
利用する
以下の2種類がある
基本プロファイラ
 主な用途：プログラム全体で、最も時間のかかって
いる関数を同定する
詳細プロファイラ
 主な用途：最も時間のかかっている関数内の指定
部分において、メモリアクセス効率、キャッシュヒット率、
スレッド実行効率、ＭＰＩ通信頻度解析、を行う
20
2023年度
計算科学技術特論A


# Page. 21

![Page Image](https://bcdn.docswell.com/page/Y7656K3PEV.jpg)

性能プロファイラの種類の詳細

基本プロファイラ





コマンド例：fipp –C
表示コマンド：fipppx、GUI（ＷＥＢ経由）
ユーザプログラムに対し一定間隔(デフォルト時100 ミリ秒間隔)毎に割り込み
をかけ情報を収集する。
収集した情報を基に、コスト情報等の分析結果を表示。
詳細プロファイラ




コマンド例：fapp –C
表示コマンド：GUI(WEB経由)
ユーザプログラムの中に測定範囲を設定し、測定範囲のハードウェア
カウンタの値を収集。
収集した情報を基に、MFLOPS、MIPS、各種命令比率、キャッシュミス等の
詳細な分析結果を表示。
21
2023年度
計算科学技術特論A


# Page. 22

![Page Image](https://bcdn.docswell.com/page/G75Y6PLP74.jpg)

基本プロファイラ利用例




調べるべきプログラムのあるディレクトリに
Profディレクトリを作成
$ mkdir Prof
wａ２(対象の実行可能ファイル) の wa2-pure.bash（スクリプト
ファイル） 中に以下を記載
fipp -C -d Prof mpirun ./wa2
実行する
$ pjsub wa2-pure.bash
テキストプロファイラを起動
$ fipppx –A -d Prof
22
2023年度
計算科学技術特論A


# Page. 23

![Page Image](https://bcdn.docswell.com/page/9J2GM6LZJR.jpg)

基本プロファイラ出力例（1/2）
-----------------------------------------------------------------------------------------Fujitsu Instant Profiler Version 1.2.0
Measured time
: Thu Apr 19 09:32:18 2012
CPU frequency
: Process
Type of program
: MPI
0-
127 1848 (MHz)
Average at sampling interval : 100.0 (ms)
Measured range
: All ranges
Virtual coordinate
: (12, 0, 0)
-----------------------------------------------------------------------------------------____________________________________________________________________________________
Time statistics
Elapsed(s)
User(s)
System(s)
--------------------------------------------2.1684
53.9800
87.0800 Application
--------------------------------------------2.1684
0.5100
0.6400 Process
11
2.1588
0.4600
0.6800 Process
88
2.1580
0.5000
0.6400 Process
99
2.1568
0.6600
1.4200 Process
111
各MPIプロセスの
経過時間、ユーザ時間、システム時間
…
23
2023年度
計算科学技術特論A


# Page. 24

![Page Image](https://bcdn.docswell.com/page/DEYVX9QNEM.jpg)

基本プロファイラ出力例（2/2）
__________________________________________________________________________________________
Procedures profile
各関数の実行時間が、
全体時間に占める割合
*************************************************************************************
Application - procedures
*************************************************************************************
Cost
%
Mpi
%
Start
End
具体的な箇所と、
ソースコード上の
行数の情報
---------------------------------------------------------------------475 100.0000
312
65.6842
--
-- Application
---------------------------------------------------------------------312
65.6842
312 100.0000
1
45 MAIN__
82
17.2632
0
0.0000
--
-- __GI___sched_yield
80
16.8421
0
0.0000
--
-- __libc_poll
1
0.2105
0
0.0000
--
-- __pthread_mutex_unlock_usercnt
*************************************************************************************
Process
11 - procedures
*************************************************************************************
Cost
%
Mpi
%
Start
End
---------------------------------------------------------------------5 100.0000
4
80.0000
--
-- Process
11
---------------------------------------------------------------------4
80.0000
4 100.0000
1
20.0000
0
0.0000
1
--
45 MAIN__
-- __GI___sched_yeld
….
24
2023年度
計算科学技術特論A


# Page. 25

![Page Image](https://bcdn.docswell.com/page/VJNM1LKV78.jpg)

詳細プロファイラ利用例


測定したい対象に、以下のコマンドを挿入
Ｆｏｒｔｒａｎ言語の場合





ヘッダファイル：なし
測定開始 手続き名： call fapp_start(name, number, level)
測定終了 手続き名： call fapp_stop(name, number, level)
利用例： call fapp_start(“region１”,１,１)
Ｃ/C++言語の場合




25
ヘッダファイル： fj_tool/fjcoll.h
測定開始 関数名：void fapp_start(const char *name, int number, int level)
測定終了 関数名： void fapp_stop(const char *name, int number, int level)
利用例： fapp_start(“region１”,１,１);
2023年度
計算科学技術特論A


# Page. 26

![Page Image](https://bcdn.docswell.com/page/YE9N54VVE3.jpg)

詳細プロファイラ利用例
/Wa2 に Profディレクトリを作成
$ mkdir Prof
 Ｗａ２のwa2-pure.bash中に以下を記載
(キャッシュ情報取得時)

fapp -C -d Prof –L 1 –Ihwm –Hevent=Cache mpirun ./wa2

実行する
$ pjsub wa2-pure.bash
26
2023年度
計算科学技術特論A


# Page. 27

![Page Image](https://bcdn.docswell.com/page/GE8KYQP47D.jpg)

CPU解析レポート（エクセル形式）






性能プロファイルは見にくい
性能プロファイルデータ（マシン語命令の種類や、実行
時間に占める割合など）を、Ｅｘｃｅｌで可視化してくれる
ツール
コマンド例：fapp –c –Hevent=pa1 ./a.out
単体レポート： １回測定
標準レポート：１１回測定
詳細レポート：１７回測定
27
2023年度
計算科学技術特論A


# Page. 28

![Page Image](https://bcdn.docswell.com/page/LELZ4X1VJR.jpg)

CPU解析レポート（エクセル形式）
 手順
1.
2.
3.
4.
5.
6.
28
対象箇所（ループ）を、専用のAPIで指定する
プロファイルを入れるフォルダを＜測定数分＞か所を
つくる
プロファイルのためのコマンドで＜測定数分＞回実
行する
エクセル形式に変換する
4のエクセル形式を手元のパソコンに持ってくる
5のファイルを、指定のエクセルと同一のフォルダに
入れてから、指定のエクセルを開く
2023年度
計算科学技術特論A


# Page. 29

![Page Image](https://bcdn.docswell.com/page/4JMVKLWNJW.jpg)

CPU解析レポートのための指示ＡＰＩ

以下のAPIで、対象となるループを挟む（Fortranの場合）
call fapp_start (“region”, １)
&lt;対象となるループ&gt;
call fapp_stop (“region”, １)



29
詳細プロファイラの指定APIと同じです
“region”は、対象となる場所の名前なので、任意の名前を付
けることが可能（後で、専用エクセルを開くときに使う）
“１”は、レベルの指定で、数字を書く
 -L オプションで指定したレベル以上を測定
2023年度
計算科学技術特論A


# Page. 30

![Page Image](https://bcdn.docswell.com/page/PJR2XKY4E9.jpg)

表示例
ソース：FUJITSU Software Technical Computing Suite V4.0L20
Development Studioプロファイラ使用手引書, J2UL-2483-02Z0(00), 2020年3月
30
2023年度
計算科学技術特論A


# Page. 31

![Page Image](https://bcdn.docswell.com/page/PEXD6LGZ7X.jpg)

そのほかの最適化技法
ループ分割、ループ消滅とスレッド並列化、等
31
2023年度
計算科学技術特論A


# Page. 32

![Page Image](https://bcdn.docswell.com/page/3EKX1L4V7D.jpg)

ループ分割・ループ消滅の事例
32
2023年度
計算科学技術特論A


# Page. 33

![Page Image](https://bcdn.docswell.com/page/L735Q3DQE5.jpg)

ループ分割とループ消滅の実例（その１）
 Seism3D：



東京大学古村教授が開発した地震波のシミュレーション
プログラム（における、ベンチマークプログラム）
 東京大学情報基盤センターで開発中の
数値計算ミドルウェアppOpen-HPCにおける
ppOpen-APPL/FDMとして開発
有限差分法（Finite Differential Method (FDM)）
３次元シミュレーション
 ３次元配列が確保される
データ型: 単精度 (real*4)
33
2023年度
計算科学技術特論A


# Page. 34

![Page Image](https://bcdn.docswell.com/page/87DD94PW7G.jpg)

ループ分割とループ消滅の実例（その２）

作業領域が多い実装のため
最大問題サイズ: NX=256, NY=256, NZ=128（32GBメモリ）
 たった 32.1MB分しか問題空間として確保できない
⇒ほとんどのデータは、キャッシュに載ってしまう
 時間ステップ数が大きい場合、総合時間を考慮すると１ステップ
当たりの実行時間を減らすしかない ⇒問題サイズを小さくする
近年のマルチコア計算機の傾向
 L3キャッシュ（Last Level Cache, LLC）が大きくなってきている
 Xeon E5-2600 v3（Haswell-EP）、LLC: 45MB [L3/socket]
 FX100、Sparc64 XI-fx、LLC: 24 MB [L2/node]
 次世代デバイス（Intel 3D Xpoint、ストレージクラスメモリ）など
ではTB級（予想）
⇒問題空間の配列データが小さい時、キャッシュ上にデータが
のりやすくなってきている
34


2023年度
計算科学技術特論A


# Page. 35

![Page Image](https://bcdn.docswell.com/page/VJPXZMNX78.jpg)

主要カーネル（第1位）: 全体の9.8%
subroutine ppohFDM_update_stress (ファイル名：m_stress.f90)
do k = NZ00, NZ01
do j = NY00, NY01
do i = NX00, NX01
RL1 = LAM (I,J,K)
RM1 = RIG (I,J,K)
RM2 = RM1 + RM1
RLRM2 = RL1+RM2
DXVX1 = DXVX(I,J,K)
DYVY1 = DYVY(I,J,K)
DZVZ1 = DZVZ(I,J,K)
D3V3 = DXVX1 + DYVY1 + DZVZ1
DXVYDYVX1 = DXVY(I,J,K)+DYVX(I,J,K)
DXVZDZVX1 = DXVZ(I,J,K)+DZVX(I,J,K)
DYVZDZVY1 = DYVZ(I,J,K)+DZVY(I,J,K)
SXX (I,J,K) = SXX (I,J,K) + (RLRM2*(D3V3)-RM2*(DZVZ1+DYVY1) ) * DT
SYY (I,J,K) = SYY (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DZVZ1) ) * DT
SZZ (I,J,K) = SZZ (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DYVY1) ) * DT
SXY (I,J,K) = SXY (I,J,K) + RM1 * DXVYDYVX1 * DT
SXZ (I,J,K) = SXZ (I,J,K) + RM1 * DXVZDZVX1 * DT
SYZ (I,J,K) = SYZ (I,J,K) + RM1 * DYVZDZVY1 * DT
end do
end do
35
2023年度 計算科学技術特論A
end do


# Page. 36

![Page Image](https://bcdn.docswell.com/page/2EVLW9K3JQ.jpg)

主要カーネル（第2位）: 全体の6.8%
subroutine ppohFDM_update_stress_sponge (ファイル名：m_stress.f90)
do k = NZ00, NZ01
gg_z = gz(k)
do j = NY00, NY01
gg_y = gy(j)
gg_yz = gg_y * gg_z
do i = NX00, NX01
gg_x = gx(i)
gg_xyz = gg_x * gg_yz
SXX(I,J,K) = SXX(I,J,K) * gg_xyz
SYY(I,J,K) = SYY(I,J,K) * gg_xyz
SZZ(I,J,K) = SZZ(I,J,K) * gg_xyz
SXY(I,J,K) = SXY(I,J,K) * gg_xyz
SXZ(I,J,K) = SXZ(I,J,K) * gg_xyz
SYZ(I,J,K) = SYZ(I,J,K) * gg_xyz
end do
end do
end do
36
2023年度
計算科学技術特論A


# Page. 37

![Page Image](https://bcdn.docswell.com/page/57GXDZQYEL.jpg)

主要カーネル（第3位）: 全体の6.2%
subroutine ppohFDM_update_vel (ファイル名：m_velocity.f90)
do k = NZ00, NZ01
do j = NY00, NY01
do i = NX00, NX01
! Effective Density
ROX = 2.0_PN/( DEN(I,J,K) + DEN(I+1,J,K) )
ROY = 2.0_PN/( DEN(I,J,K) + DEN(I,J+1,K) )
ROZ = 2.0_PN/( DEN(I,J,K) + DEN(I,J,K+1) )
VX(I,J,K) = VX(I,J,K) + ( DXSXX(I,J,K)+DYSXY(I,J,K)+DZSXZ(I,J,K) )*ROX*DT
VY(I,J,K) = VY(I,J,K) + ( DXSXY(I,J,K)+DYSYY(I,J,K)+DZSYZ(I,J,K) )*ROY*DT
VZ(I,J,K) = VZ(I,J,K) + ( DXSXZ(I,J,K)+DYSYZ(I,J,K)+DZSZZ(I,J,K) )*ROZ*DT
end do
end do
end do
37
2023年度
計算科学技術特論A


# Page. 38

![Page Image](https://bcdn.docswell.com/page/4EQ81LK6JP.jpg)

主要カーネル（第4位）: 全体の5.8%
subroutine ppohFDM_pdiffy3_p4 (ファイル名：m_pfd3d.f90)
R40 = C40/DY
R41 = C41/DY
do K = 1, NZ
do I = 1, NX
do J = 1, NY
DYV (I,J,K) = (V(I,J+1,K)-V(I,J,K) )*R40 - (V(I,J+2,K)-V(I,J-1,K))*R41
end do
end do
end do
38
2023年度
計算科学技術特論A


# Page. 39

![Page Image](https://bcdn.docswell.com/page/KJ42RDPME1.jpg)

カーネルループの構造

以下の３重ループを検討する
（ppOpen-APPL/FDMの第１位ループと同等）
DO K = １, NZ
DO J = １, NY
DO I = １，NX
RL = LAM (I,J,K)
RM = RIG (I,J,K)
RM2 = RM + RM
RMAXY = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I,J+1,K+1))
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
END DO
END DO
END DO
39
2023年度
計算科学技術特論A


# Page. 40

![Page Image](https://bcdn.docswell.com/page/LE1MQZ9Y7G.jpg)

ここでのコード最適化の方針（その１）
「富岳」ではレジスタ数
が少なく、スピルコード
が出やすい
→ループ分割で回避
(Loop Splitting)
 スピルコード
（レジスタから追い出されるデータがある
コード）を防ぐ目的で行う。
 レジスタを最大限に使うプログラムで、
メモリからのデータ読み出しを削減し、
高速化する。
 ループ分割
40
2023年度
計算科学技術特論A


# Page. 41

![Page Image](https://bcdn.docswell.com/page/GEW3M9Q1E2.jpg)

ここでのコード最適化の方針（その２）

ループ消滅（Loop Collapse)
 対象は３重ループ → 以下の２つの方針がある
 １次元ループ化
 スレッド並列実行のため、最外側のループ長を増加させる
目的で行う
 ベクトル計算機用のコンパイラで行われることが多い
 メニーコア計算機でも状況により効果が見込まれる
 ２次元ループ化
 スレッド並列実行のため、最外側のループ長を増加させる
目的で行う
 コンパイラによる最内ループのプリフェッチ処理を増進
 近年のメニーコア計算機でもっとも有望と思われる方法
41
2023年度
計算科学技術特論A


# Page. 42

![Page Image](https://bcdn.docswell.com/page/47ZVQ9WXE3.jpg)

ループ分割の例 – 分割点

例:以下の箇所でループ分割する例
DO K = １, NZ
DO J = １, NY
DO I = １，NX
RL(I) = LAM (I,J,K)
RM(I) = RIG (I,J,K)
RM2(I) = RM(I) + RM(I)
RMAXY(I) = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ(I) = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I+1,J,K+1))
RMAYZ(I) = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I,J+1,K+1))
RLTHETA(I) = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL(I)
QG(I) = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
END DO
ループ分割点
DO I = １, NX
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA(I) + RM2(I)*DXVX(I,J,K))*DT )*QG(I)
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA(I) + RM2(I)*DYVY(I,J,K))*DT )*QG(I)
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA(I) + RM2(I)*DZVZ(I,J,K))*DT )*QG(I)
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY(I)*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG(I)
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ(I)*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG(I)
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ(I)*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG(I)
END DO
END DO
42
2023年度 計算科学技術特論A
END
DO


# Page. 43

![Page Image](https://bcdn.docswell.com/page/YJ656KGP7V.jpg)

ループ消滅 – １重ループ化

例）
利点：ループ長が増える
NZ → NZ*NY*NX
DO KK = １, NZ * NY * NX
K = (KK-１)/(NY*NX) + １
J = mod((KK-１)/NX,NY) + １
I = mod(KK-１,NX) + １
RL = LAM (I,J,K)
RM = RIG (I,J,K)
RM2 = RM + RM
RMAXY = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I,J+1,K+1))
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
END DO
43
2023年度
計算科学技術特論A


# Page. 44

![Page Image](https://bcdn.docswell.com/page/GJ5Y6PVPJ4.jpg)

ループ消滅 – ２重ループ化

例）
利点：ループ長が増える
NZ → NZ*NY
DO KK = １, NZ * NY
K = (KK-１)/NY + １
J = mod(KK-１,NY) + １
DO I = １, NX
このI-ループは連続：
RL = LAM (I,J,K)
コンパイラによるプリフェッチコード生成が可能
RM = RIG (I,J,K)
RM2 = RM + RM
RMAXY = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I,J+1,K+1))
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
ENDDO
END DO
44
2023年度
計算科学技術特論A


# Page. 45

![Page Image](https://bcdn.docswell.com/page/LE35Q3D6J5.jpg)

さらなる改良：定義－参照距離の変更
DO K = １, NZ
DO J = １, NY
DO I = １, NX
RL = LAM (I,J,K)
RM = RIG (I,J,K)
RM2 = RM + RM
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K)+ (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K)+ (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
RMAXY = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I+1,J,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(1.0/RIG(I,J,K) + 1.0/RIG(I,J+1,K) + 1.0/RIG(I,J,K+1) + 1.0/RIG(I,J+1,K+1))
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
END DO
END DO
END DO
T2K（AMD Opteron）で、約50%の速度向上
45
2023年度
計算科学技術特論A


# Page. 46

![Page Image](https://bcdn.docswell.com/page/8EDD94PMEG.jpg)

修正コード + I-ループ分割の例
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
RL = LAM (I,J,K)
RM = RIG (I,J,K)
RM2 = RM + RM
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
ENDDO
DO I = 1, NX
STMP1 = 1.0/RIG(I,J,K)
STMP2 = 1.0/RIG(I+1,J,K)
STMP4 = 1.0/RIG(I,J,K+1)
STMP3 = STMP1 + STMP2
RMAXY = 4.0/(STMP3 + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I,J+1,K+1))
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
END DO
END DO
END
46 DO
2023年度 計算科学技術特論A
ループ分割すると、
QGの再計算が必要になる
通常のコンパイラでは
ユーザの判断が必要
なので、できない


# Page. 47

![Page Image](https://bcdn.docswell.com/page/V7PXZMNQE8.jpg)

修正コード + K-ループの分割の例
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
RL = LAM (I,J,K)
RM = RIG (I,J,K)
RM2 = RM + RM
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
ENDDO; ENDDO; ENDDO
完全に別の３重ループに分かれる
←分かれた３重ループに対し、
コンパイラによるさらなる最適化の可能性
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
STMP1 = 1.0/RIG(I,J,K)
STMP2 = 1.0/RIG(I+1,J,K)
STMP4 = 1.0/RIG(I,J,K+1)
STMP3 = STMP1 + STMP2
RMAXY = 4.0/(STMP3 + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I,J+1,K+1))
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
47
2023年度 計算科学技術特論A
END DO; END DO; END DO;


# Page. 48

![Page Image](https://bcdn.docswell.com/page/2JVLW9KP7Q.jpg)

チューニングの可能性のあるコード例
（経験的に決めた数例について）







#1 : 基の３重ループコード（ベースライン）
#2: I-ループ分割のみ
#3: J-ループ分割のみ
#4: K-ループ分割のみ
#5: #2ループに対するループ消滅 (2重ループ化)
#6 : #1ループに対するループ消滅 (１重ループ化)
#7 : #1ループに対するループ消滅 (2重ループ化)
48
2023年度
計算科学技術特論A


# Page. 49

![Page Image](https://bcdn.docswell.com/page/5EGXDZQQJL.jpg)

ループ分割・ループ消滅の効果

東京大学情報基盤センターFX１０を利用



最外ループに対して、OpenMPが適用可能


49
１ノード、１６スレッド
Sparc64 IV-fx (1.848 GHz)
parallel do構文で並列化可能
スレッド数は、１～１６まで変更可能
2023年度
計算科学技術特論A


# Page. 50

![Page Image](https://bcdn.docswell.com/page/4JQ81LKW7P.jpg)

チューニングの結果
48.78
実行時間[秒]
25
#1
#2
#3
#5
#6
#7
#4
20
15
10
5
0
1
2
スレッド数
50
4
8
#2
16
2023年度
#3
#4
#1
計算科学技術特論A
#5
#6
#7
実装の
種類
＃４の実装が常に
最高速


# Page. 51

![Page Image](https://bcdn.docswell.com/page/K742RDP1J1.jpg)

＃１（ベースライン）に対する速度向上
SpeedUP
1.8
1.6
1.54
1.54
1.53
1.53
1
2
4
8
1.52
1.4
1.2
1
0.8
0.6
0.4
0.2
0
51
2023年度
計算科学技術特論A
16 スレッド数


# Page. 52

![Page Image](https://bcdn.docswell.com/page/LJ1MQZ95EG.jpg)

#4のK-ループの分割のコード
!$omp parallel do private(k,j,i,STMP1,STMP2,STMP3,STMP4,RL,RM,RM2,
!$omp&amp; RMAXY,RMAXZ,RMAYZ,RLTHETA,QG)
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
RL = LAM (I,J,K); RM = RIG (I,J,K); RM2 = RM + RM;
RLTHETA = (DXVX(I,J,K)+DYVY(I,J,K)+DZVZ(I,J,K))*RL
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXX (I,J,K) = ( SXX (I,J,K) + (RLTHETA + RM2*DXVX(I,J,K))*DT )*QG
SYY (I,J,K) = ( SYY (I,J,K) + (RLTHETA + RM2*DYVY(I,J,K))*DT )*QG
SZZ (I,J,K) = ( SZZ (I,J,K) + (RLTHETA + RM2*DZVZ(I,J,K))*DT )*QG
ENDDO; ENDDO; ENDDO
!$omp end parallel do
!$omp parallel do private(k,j,i,STMP1,STMP2,STMP3,STMP4,RL,RM,RM2,
!$omp&amp; RMAXY,RMAXZ,RMAYZ,RLTHETA,QG)
DO K = 1, NZ
DO J = 1, NY
DO I = 1, NX
STMP1 = 1.0/RIG(I,J,K); STMP2 = 1.0/RIG(I+1,J,K); STMP4 = 1.0/RIG(I,J,K+1);
STMP3 = STMP1 + STMP2
RMAXY = 4.0/(STMP3 + 1.0/RIG(I,J+1,K) + 1.0/RIG(I+1,J+1,K))
RMAXZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I+1,J,K+1))
RMAYZ = 4.0/(STMP3 + STMP4 + 1.0/RIG(I,J+1,K+1))
QG = ABSX(I)*ABSY(J)*ABSZ(K)*Q(I,J,K)
SXY (I,J,K) = ( SXY (I,J,K) + (RMAXY*(DXVY(I,J,K)+DYVX(I,J,K)))*DT )*QG
SXZ (I,J,K) = ( SXZ (I,J,K) + (RMAXZ*(DXVZ(I,J,K)+DZVX(I,J,K)))*DT )*QG
SYZ (I,J,K) = ( SYZ (I,J,K) + (RMAYZ*(DYVZ(I,J,K)+DZVY(I,J,K)))*DT )*QG
END DO; END DO; END DO;
52
2023年度 計算科学技術特論A
!$omp end parallel do


# Page. 53

![Page Image](https://bcdn.docswell.com/page/GJW3M9QWJ2.jpg)

マルチコアCPUではループ消滅が有効？



一般に、３次元陽解法のカーネルは以下の構造
OpenMPのスレッド並列化は最外側ループに適用
この時、並列性はK-ループ長のNZで決まる
!$omp parallel do private(…) （NZ＞スレッド数） が並列性のため必要
DO K = １, NZ
 OpenMPオーバヘッドを考えると、ノードあたりのNZ
DO J = １, NY
はスレッド数の２～３倍必要
DO I = １, NX
 例）68スレッドなら、NZは140～210以上
＜離散化手法に基づく数式＞  HTで272スレッド実行なら、NZは540～810以上
 3次元問題で上記のサイズ（ノード当たりの
ENDDO
問題サイズ）は確保できるか？
ENDDO
ENDDO
ループ長が確保できない場合
!$omp end parallel do
ループ消滅が必須
53
2023年度
計算科学技術特論A


# Page. 54

![Page Image](https://bcdn.docswell.com/page/4EZVQ9W2J3.jpg)

計算機環境
(Xeon Phi (KNC) の８ノード)
 Intel
Xeon Phi (KNC)

Xeon Phi 5110P (1.053 GHz), 60 cores

メモリ量：8 GB (GDDR5)
理論ピーク性能：1.01 TFLOPS
Xeon Phiのクラスタ（ノード当たり１ボード）
InfiniBand FDR x 2 Ports









Intel MPI





Mellanox Connect-IB
PCI-E Gen3 x16
56Gbps x 2
理論バンド幅 13.6GB/s
フルバイセクション
MPICH2、MVAPICH2ベース
4.1 Update 3 (build 048)
コンパイラ：Intel Fortran version 14.0.0.080 Build 20130728
コンパイラオプション：
-ipo20 -O3 -warn all -openmp -mcmodel=medium -shared-intel -mmic -align array64byte
KMP_AFFINITY=granularity=fine, balanced (ソケット間へスレッドを均等割り当て)
54
2023年度
計算科学技術特論A


# Page. 55

![Page Image](https://bcdn.docswell.com/page/Y7656KX6EV.jpg)

実行詳細
ppOpen-APPL/FDM ver.0.2
 ppOpen-AT ver.0.2
 時間ステップ数: 2000 steps
 ノード数: 8 ノード
 Native Mode 実行
 問題サイズ
(8 GB/ノードでの最大サイズ)



55
NX * NY * NZ = 1536 x 768 x 240 / 8ノード
NX * NY * NZ = 768 * 384 * 120 / ノード
(MPIプロセス当たりのサイズではない)
2023年度
計算科学技術特論A


# Page. 56

![Page Image](https://bcdn.docswell.com/page/G75Y6PX574.jpg)

ハイブリッドMPI/OpenMP実行の詳細

Xeon PhiにおけるMPIプロセス数とOpenMPスレッド数
 4 HT (Hyper Threading)
 PX








TY: X MPIプロセス、 Y スレッド/プロセス
P8T240 : 最少のハイブリッドMPI/OpenMP実行
（ppOpen-APPL/FDMでは 8MPIプロセスが最低でも必要なため）
P16T120
# # # # # #
# # # # # # # # # #
P32T60
P2T8 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
P64T30
1つの
# # # # # #
P128T15
# # # # # # # # # #
P4T4 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 MPIプロセス
で割り当て対象
P240T8
P480T4
P960T2以下は、 この環境ではMPIエラーで実行できなかったため
除外
56
2023年度
計算科学技術特論A


# Page. 57

![Page Image](https://bcdn.docswell.com/page/9J2GM6XGJR.jpg)

スレッド当たりのループ長 (Z-軸)
(8ノード、1536x768x240 / 8 ノード )
Loop length per thread
12
プログラム上の制約で、
それぞれのMPI/OpenMP実行で
Z-軸のMPIプロセスは異なる
6
4
0.5
57
1
2
2023年度
2
計算科学技術特論A


# Page. 58

![Page Image](https://bcdn.docswell.com/page/DEYVX98GEM.jpg)

ループ消滅等による
最大の速度向上(Xeon Phi、8ノード)
Speedup [%]
558
Speedup =
max ( オリジナルコードの実行時間 / 最速実装での実行時間 )
：すべてのハイブリッドMPI/OpenMP実行 (PXTY)において
NX*NY*NZ = 1536x768x240/ 8ノード
200
171
30
20
51
演算カーネルの種類
58
2023年度
計算科学技術特論A


# Page. 59

![Page Image](https://bcdn.docswell.com/page/VJNM1LQ878.jpg)

最速のコード (update_stress)
Xeon
Phi(KNC) (P240T8)
!$omp parallel do private
(k,j,i,RL1,RM1,RM2,RLRM2,DXVX1,DYVY1,DZVZ1,D3V3,DXVYDYVX1,DXVZDZVX1,DYVZDZV1)
DO k_j = 1 , (NZ01-NZ00+1)*(NY01-NY00+1)
k = (k_j-1)/(NY01-NY00+1) + NZ00
j = mod((k_j-1),(NY01-NY00+1)) + NY00
ループ融合により
ループ長
（＝並列性）が増加
DO i = NX00, NX01
RL1 = LAM (I,J,K); RM1 = RIG (I,J,K)
RM2 = RM1 + RM1; RLRM2 = RL1+RM2
DXVX1 = DXVX(I,J,K); DYVY1 = DYVY(I,J,K); DZVZ1 = DZVZ(I,J,K);
D3V3 = DXVX1 + DYVY1 + DZVZ1;
SXX (I,J,K) = SXX (I,J,K) + (RLRM2*(D3V3)-RM2*(DZVZ1+DYVY1) ) * DT
SYY (I,J,K) = SYY (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DZVZ1) ) * DT
SZZ (I,J,K) = SZZ (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DYVY1) ) * DT
DXVYDYVX1 = DXVY(I,J,K)+DYVX(I,J,K);
DXVZDZVX1 = DXVZ(I,J,K)+DZVX(I,J,K);
DYVZDZVY1 = DYVZ(I,J,K)+DZVY(I,J,K)
SXY (I,J,K) = SXY (I,J,K) + RM1 * DXVYDYVX1 * DT
SXZ (I,J,K) = SXZ (I,J,K) + RM1 * DXVZDZVX1 * DT
SYZ (I,J,K) = SYZ (I,J,K) + RM1 * DYVZDZVY1 * DT
END DO
END DO
59
!$omp end parallel do
2023年度
計算科学技術特論A


# Page. 60

![Page Image](https://bcdn.docswell.com/page/YE9N54DXE3.jpg)

SIMD化・ループ分割の例
（「富岳」）
60
2023年度
計算科学技術特論A


# Page. 61

![Page Image](https://bcdn.docswell.com/page/GE8KYQ397D.jpg)

SIMD化・ループ分割例（「富岳」）
フラグメント分子軌道 (FMO) 法 のプログラム
ABINIT-MP [1]
 2電子積分ルーチンのSIMD化とループ分割を
実施 （2021/12/9版）


以下の関数に実施
（ SSSS, PSSS, SPSS, SSPS, SSSP, PPSS, PSPS, PSSP,
SPPS, SPSP, SSPP, DSSS, SDSS, SSDS, SSSD ）
[1] 望月祐志, 中野達也, 坂倉耕太, 渡邊啓正, 佐藤伸哉, 奥脇弘次, 秋澤和輝, 土居英男, 大島聡史,
片桐孝洋: FMOプログラムABINIT-MPの整備状況2022, Journal of Computer Chemistry, Japan,
2022年21巻4号, pp.106-110
61
2023年度
計算科学技術特論A


# Page. 62

![Page Image](https://bcdn.docswell.com/page/LELZ4XVMJR.jpg)

立教大学 望月祐志教授提供
SIMD化した積分ルーチンの例（sssp)
subroutine sub_sssp(zetam,pm,dkabm,etam,qm,dkcdm, &amp;
ma,mb,mc,md,ngij,ngkl,a,b,c,d,sint,tv)
!
!
Nov.05,&#039;02
!
T.NAKANO &amp; Y. ABE
!
use constant
use auxiliary_integral_table
use integral_parameter
implicit none
real(8),intent(in)::zetam(*),pm(3,*),dkabm(*), &amp;
etam(*),qm(3,*),dkcdm(*)
integer,intent(in)::ma,mb,mc,md,ngij,ngkl
real(8),intent(in)::a(3),b(3),c(3),d(3),tv
real(8),intent(out)::sint(*)
!---------------------------------------------integer npq,nrs,ix
real(8) p(3),q(3),qd(3),pq(3),wq(3),f(0:max_m), &amp;
dkab,zeta,dkcd,eta,ze,rz,re,rho,a0,tt
integer ts,i,j,k,l,m
real(8) delta,t_inv
real(8) ssss(0:1),f0,f1,qd1,qd2,qd3,wq1,wq2,wq3
do npq=1,ngij
if (abs(dkabm(npq)) &gt; tv) then
do nrs=1,ngkl
if (abs(dkabm(npq)*dkcdm(nrs)) &gt; tv) then
ze = 1.0_8/(zetam(npq)+etam(nrs))
a0 = dkabm(npq)*dkcdm(nrs)*sqrt(ze)
rz = etam(nrs)*ze
re = zetam(npq)*ze
rho = zetam(npq)*rz
!
!
do i=1,3
qd(i) = qm(i,nrs)-d(i)
pq(i) = qm(i,nrs)-pm(i,npq)
wq(i) =-re*pq(i)
end do
qd1
qd2
qd3
wq1
wq2
wq3
= qm(1,nrs)-d(1)
= qm(2,nrs)-d(2)
= qm(3,nrs)-d(3)
=-re*pq(1)
=-re*pq(2)
=-re*pq(3)
sint(1:3) = 0.0_8
2023/5/13
!ocl
!ocl
!ocl
!ocl
!ocl
!ocl
eval
fp_relaxed
fp_contract
noswp
eval_concurrent
SIMD
以下、次頁
62


# Page. 63

![Page Image](https://bcdn.docswell.com/page/4JMVKLG6JW.jpg)

立教大学 望月祐志教授提供
SIMD化した積分ルーチンの例（続き）
tt = (pq(1)*pq(1)+pq(2)*pq(2)+pq(3)*pq(3))*rho
if (tt &lt;= 38.0_8) then ! Tf = 2*m+36 (for m=1)
ts
= 0.5_8+tt*fmt_inv_step_size
delta = ts*fmt_step_size-tt
!
!
!
!
!
!
!
!
!
!
2023/5/13
f(0) = ((fmt_table(3,ts)*inv6*delta
+
fmt_table(2,ts)*inv2)*delta
+
fmt_table(1,ts))*delta
+
fmt_table(0,ts)
f(1) = ((fmt_table(4,ts)*inv6*delta
+
fmt_table(3,ts)*inv2)*delta
+
fmt_table(2,ts))*delta
+
fmt_table(1,ts)
f0
= ((fmt_table(3,ts)*inv6*delta
+
fmt_table(2,ts)*inv2)*delta
+
fmt_table(1,ts))*delta
+
fmt_table(0,ts)
f1
= ((fmt_table(4,ts)*inv6*delta
+
fmt_table(3,ts)*inv2)*delta
+
fmt_table(2,ts))*delta
+
fmt_table(1,ts)
else
t_inv = inv2/tt
f(0) = sqrt(pi_over2*t_inv)
f(1) = t_inv*f(0)
f0
= sqrt(pi_over2*t_inv)
f1
= t_inv*f0
end if
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
&amp;
!---------------------------------------------! ERI code generator Ver.20020228
! 2002/02/28
! T. Nakano
!
! (sssp)
!
!
ssss(0:1)=f(0:1)*a0
ssss(0)=f0*a0
ssss(1)=f1*a0
!
do l=1, 3
!
sint(l) = sint(l)+qd(l)*ssss(0)+wq(l)*ssss(1)
!
end do
sint(1) = sint(1)+qd1*ssss(0)+wq1*ssss(1)
sint(2) = sint(2)+qd2*ssss(0)+wq2*ssss(1)
sint(3) = sint(3)+qd3*ssss(0)+wq3*ssss(1)
!---------------------------------------------end if
end do
end if
end do
end subroutine sub_sssp
63


# Page. 64

![Page Image](https://bcdn.docswell.com/page/PJR2XKZNE9.jpg)

立教大学 望月祐志教授提供
改良ssppルーチン#1
subroutine sub_sspp(zetam,pm,dkabm,etam,qm,dkcdm, &amp;
ma,mb,mc,md,ngij,ngkl,a,b,c,d,sint,tv)
!
!
Nov.06,&#039;02
!
T.NAKANO &amp; Y. ABE
!
use constant
use auxiliary_integral_table
use integral_parameter
implicit none
real(8),intent(in)::zetam(*),pm(3,*),dkabm(*), &amp;
etam(*),qm(3,*),dkcdm(*)
integer,intent(in)::ma,mb,mc,md,ngij,ngkl
real(8),intent(in)::a(3),b(3),c(3),d(3),tv
real(8),intent(out)::sint(*)
!-----------------------------------------------integer npq,nrs,i,k,m,ix
real(8) :: pq(3),ze,rz,re,rho,tt
integer ts
real(8) delta,t_inv,zssss
real(8)::f0,f1,f2
real(8) ssss(0:2),ssps(1:3,0:1),sspp(1:3,1:3,0:0)
real(8) :: xqc(ngij*ngkl,3),xqd(ngij*ngkl,3),xwq(ngij*ngkl,3)
real(8) :: xtt(ngij*ngkl),xre(ngij*ngkl),xeta2(ngij*ngkl),xa0(ngij*ngkl)
integer :: npqrs
2023/5/13
64


# Page. 65

![Page Image](https://bcdn.docswell.com/page/PEXD6L987X.jpg)

立教大学 望月祐志教授提供
改良ssppルーチン#2
ix = 0
!ocl eval
!ocl fp_relaxed
!ocl fp_contract
!ocl noswp
!ocl eval_concurrent
!ocl SIMD
do npq=1,ngij
if (abs(dkabm(npq)) &lt;= tv) cycle
do nrs=1,ngkl
if (abs(dkabm(npq)*dkcdm(nrs)) &lt;= tv) cycle
ix = ix + 1
ze
= 1.0_8/(zetam(npq)+etam(nrs))
rz
= etam(nrs)*ze
re
= zetam(npq)*ze
rho = zetam(npq)*rz
do i=1,3
pq(i) = qm(i,nrs)-pm(i,npq)
xqc(ix,i) = qm(i,nrs)-c(i)
xqd(ix,i) = qm(i,nrs)-d(i)
xwq(ix,i) = -re*pq(i)
end do
xtt(ix) =
xre(ix) =
xa0(ix)
xeta2(ix)
2023/5/13
enddo
enddo
(pq(1)*pq(1)+pq(2)*pq(2)+pq(3)*pq(3))*rho
re
= dkabm(npq)*dkcdm(nrs)*sqrt(ze)
= 0.5_8/etam(nrs)
65


# Page. 66

![Page Image](https://bcdn.docswell.com/page/3EKX1LPL7D.jpg)

立教大学 望月祐志教授提供
Ver. 2 Rev. 4での速度向上の例
HIV-1 protease / FMO-MP2/6-31G* / Benchmark 100 nodes @ Fugaku
Ver. 1 Rev. 22
Ver. 2 Rev. 4
=====================
## TIME PROFILE
=====================
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
time:
time:
time:
time:
time:
time:
time:
time:
=====================
## TIME PROFILE
=====================
Monomer SCF
Monomer MP2
Monomer (Total)
Dimer ES
Dimer SCF
Dimer MP2
Dimer (Total)
FMO (Total)
## Time profile
452.4
17.4
472.7
99.7
278.6
269.1
695.3
1168.0
seconds
seconds
seconds
seconds
seconds
seconds
seconds
seconds
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
Elapsed
time:
time:
time:
time:
time:
time:
time:
time:
Monomer SCF
Monomer MP2
Monomer (Total)
Dimer ES
Dimer SCF
Dimer MP2
Dimer (Total)
FMO (Total)
=
=
=
=
=
=
=
=
354.7
16.0
373.4
109.5
221.7
242.4
673.4
1046.8
seconds
seconds
seconds
seconds
seconds
seconds
seconds
seconds
## Time profile
Number of cores (total)
=
Number of cores (fragment) =
THREADS (FRAGMENT)
Total time =
=
=
=
=
=
=
=
=
=
200
1
24
1172.8 seconds
Number of cores (total)
=
Number of cores (fragment) =
THREADS (FRAGMENT)
Total time =
=
200
1
24
1050.2 seconds
「富岳」-2021年9月時点
2023/5/13
・ Ver. 2 Rev. 4はA64FX向け積分SIMD化、「不要配列」の整理などを反映済み
・ より大型の系ではMP2ジョブで2-5割程度の速度向上
・ cc-pVDZの方が短縮長が長いために加速効果が出やすい （他系でも評価）
66


# Page. 67

![Page Image](https://bcdn.docswell.com/page/L735Q3X6E5.jpg)

立教大学 望月祐志教授提供
改良ssppルーチン#3
sint(1:9)=0.0_8
!ocl eval
!ocl fp_relaxed
!ocl fp_contract
!ocl noswp
!ocl eval_concurrent
!ocl SIMD
do npqrs = 1,ix
tt = xtt(npqrs)
if (tt &lt;= 40.0_8) then ! Tf=2*m+36
ts=0.5_8+tt*fmt_inv_step_size
delta=ts*fmt_step_size-tt
f0=((fmt_table(03,ts)*inv6*delta &amp;
+fmt_table(02,ts)*inv2)*delta &amp;
+fmt_table(01,ts))*delta &amp;
+fmt_table(00,ts)
f1=((fmt_table(04,ts)*inv6*delta &amp;
+fmt_table(03,ts)*inv2)*delta &amp;
+fmt_table(02,ts))*delta &amp;
+fmt_table(01,ts)
f2=((fmt_table(05,ts)*inv6*delta &amp;
+fmt_table(04,ts)*inv2)*delta &amp;
+fmt_table(03,ts))*delta &amp;
+fmt_table(02,ts)
else
t_inv=inv2/tt
f0=sqrt(pi_over2*t_inv)
f1=t_inv*f0
f2=t_inv*3.0_8*f1
end if
2023/5/13
ssss(0)=f0*xa0(npqrs)
ssss(1)=f1*xa0(npqrs)
ssss(2)=f2*xa0(npqrs)
do m=0, 1
ssps(1,m)=xqc(npqrs,1)*ssss(m)+xwq(npqrs,1)*ssss(m+1)
ssps(2,m)=xqc(npqrs,2)*ssss(m)+xwq(npqrs,2)*ssss(m+1)
ssps(3,m)=xqc(npqrs,3)*ssss(m)+xwq(npqrs,3)*ssss(m+1)
end do
do k=1, 3
sspp(k,1,0)=xqd(npqrs,1)*ssps(k,0)+xwq(npqrs,1)*ssps(k,1)
sspp(k,2,0)=xqd(npqrs,2)*ssps(k,0)+xwq(npqrs,2)*ssps(k,1)
sspp(k,3,0)=xqd(npqrs,3)*ssps(k,0)+xwq(npqrs,3)*ssps(k,1)
end do
zssss=xeta2(npqrs)*(ssss(0)-xre(npqrs)*ssss(1))
sspp(1,1,0)=sspp(1,1,0)+zssss
sspp(2,2,0)=sspp(2,2,0)+zssss
sspp(3,3,0)=sspp(3,3,0)+zssss
sint(1)
sint(2)
sint(3)
sint(4)
sint(5)
sint(6)
sint(7)
sint(8)
sint(9)
=
=
=
=
=
=
=
=
=
sint(1)+sspp(1,1,0)
sint(2)+sspp(1,2,0)
sint(3)+sspp(1,3,0)
sint(4)+sspp(2,1,0)
sint(5)+sspp(2,2,0)
sint(6)+sspp(2,3,0)
sint(7)+sspp(3,1,0)
sint(8)+sspp(3,2,0)
sint(9)+sspp(3,3,0)
enddo
end subroutine sub_sspp
67


# Page. 68

![Page Image](https://bcdn.docswell.com/page/87DD942M7G.jpg)

立教大学 望月祐志教授提供
積分ルーチンのループ分割例
Sub_ssssの例
do npq=1,ngij
if (abs(dkabm(npq)) &gt; tv) then
do nrs=1,ngkl
if (abs(dkabm(npq)*dkcdm(nrs)) &gt; tv) then
ze = 1.0_8/(zetam(npq)+etam(nrs))
a0 = dkabm(npq)*dkcdm(nrs)*sqrt(ze)
rz = etam(nrs)*ze
rho = zetam(npq)*rz
tt = ((qm(1,nrs)-pm(1,npq))*(qm(1,nrs)-pm(1,npq)) &amp;
+(qm(2,nrs)-pm(2,npq))*(qm(2,nrs)-pm(2,npq)) &amp;
+(qm(3,nrs)-pm(3,npq))*(qm(3,nrs)-pm(3,npq)))*rho
if (tt &lt;= 36.0_8) then ! Tf = 2*m+36 (for the case of m=0)
ts = 0.5_8+tt*fmt_inv_step_size
delta = ts*fmt_step_size-tt
ssss(0) = (((fmt_table(3,ts)*inv6*delta &amp;
+fmt_table(2,ts)*inv2)*delta &amp;
+fmt_table(1,ts))*delta &amp;
+fmt_table(0,ts))*a0
else
ssss(0) = sqrt(pi_over4/tt)*a0
end if
sint(1) = sint(1)+ssss(0)
end if
end do
end if
end do
do npq=1,ngij
if (abs(dkabm(npq)) &lt;= tv) cycle
do nrs=1,ngkl
if (abs(dkabm(npq)*dkcdm(nrs)) &lt;= tv) cycle
ix = ix + 1
ze = 1.0_8/(zetam(npq)+etam(nrs))
xa0(ix) = dkabm(npq)*dkcdm(nrs)*sqrt(ze)
rz = etam(nrs)*ze
rho = zetam(npq)*rz
xtt(ix) = ((qm(1,nrs)-pm(1,npq))*(qm(1,nrs)-pm(1,npq)) &amp;
+(qm(2,nrs)-pm(2,npq))*(qm(2,nrs)-pm(2,npq)) &amp;
+(qm(3,nrs)-pm(3,npq))*(qm(3,nrs)-pm(3,npq)))*rho
enddo
enddo
sint(1) = 0.0_8
do npqrs=1,ix
tt = xtt(npqrs)
if (tt &lt;= 36.0_8) then ! Tf = 2*m+36 (for the case of m=0)
ts = 0.5_8+tt*fmt_inv_step_size
delta = ts*fmt_step_size-tt
ssss(0) = (((fmt_table(3,ts)*inv6*delta &amp;
+fmt_table(2,ts)*inv2)*delta &amp;
+fmt_table(1,ts))*delta &amp;
+fmt_table(0,ts))*xa0(npqrs)
else
ssss(0) = sqrt(pi_over4/tt)*xa0(npqrs)
end if
sint(1) = sint(1)+ssss(0)
enddo
 積分タイプにも拠るが20%～50%の加速
 ループ分割はレジスタスピルの低減に効果？
手動にてSIMD化とループ分割を併用する修正を行って40%の高速化（約1.6倍）
2023/5/13
68


# Page. 69

![Page Image](https://bcdn.docswell.com/page/VJPXZM1Q78.jpg)

そのほかの最適化技法
ハイブリッドMPI/OpenMPの実行形態
69
2023年度
計算科学技術特論A


# Page. 70

![Page Image](https://bcdn.docswell.com/page/2EVLW9YPJQ.jpg)

名古屋大学FX100 (2020年3月末退役済み)
The Fujitsu PRIMEHPC FX100
Contents
Whole
System
Specifications
Total Performance
3.2 PFLOPS
Total Memory Amounts
90 TiB
Total #nodes
2,880
Inter Connection
The TOFU2
(6 Dimension
Mesh / Torus)
Local File System Amounts
6.0 PB
2880 Nodes (92,160 Cores)
Contents
Node
Processor
Specifications
Theoretical Peak Performance
1 TFLOPS (double precision)
#Processors (#Cores)
32 + 2 (assistant cores)
Main Memory Amounts
32 GB
Processor Name
SPARC64 XI-fx
Frequency
2.2 GHz
Theoretical Peak Performance (Core)
31.25 GFLOPS
70


# Page. 71

![Page Image](https://bcdn.docswell.com/page/57GXDZ9QEL.jpg)

スーパーコンピュータ「不老」
Type I サブシステム(FX1000)
機種名
FUJITSU PRIMEHPC FX1000
計算ノード
CPU
A64FX (Armv8.2-A + SVE), 48コア+2アシスタント
コア( I/O兼計算ノードは48コア+ 4アシスタント
コア ), 2.2GHz, 1ソケット
メインメモリ
HBM2, 32GiB
倍精度 3.3792 TFLOPS, 単精度 6.7584 TFLOPS,
理論演算性能
半精度 13.5168 TFLOPS
1,024 GB/s (1CMG=12コアあたり256 GB/s,
メモリバンド幅
1CPU=4CMG)
2,304ノード, 110,592コア
ノード数、総コア数
(+4,800アシスタントコア)
7.782 PFLOPS
総理論演算性能
72 TiB
総メモリ容量
TofuインターコネクトD
各ノードは周囲の隣接ノードへ同時に
ノード間インターコネクト 合計 40.8 GB/s × 双方向 で通信可能（1リンク
当たり 6.8 GB/s × 双方向, 6リンク同時通信
可能）
ユーザ用
なし
ローカルストレージ
冷却方式
水冷
2023年度 計算科学技術特論A
ノード内構成
71


# Page. 72

![Page Image](https://bcdn.docswell.com/page/4EQ81LWWJP.jpg)

実行詳細
•
•
•
•
ppOpen-APPL/FDM ver.1.0 （差分法コード）
時間ステップ数: 2000 steps
ノード数: 8 node
問題サイズ (8 GB/node)
– NX * NY * NZ = 512 x 512 x 512 / 8 Node
– NX * NY * NZ = 256 * 256 * 256 / node (!= MPI プロセス)
• 対象 MPIプロセス数 と スレッド数
– PXTY: X MPIプロセス、 Y スレッド／プロセス
– P8T32 : 最小の MPI-OpenMP実行条件（ppOpen-APPL/FDM）：
最低でも 8 MPIプロセスが要るため
– P16T16
– P32T8
– P64T4
– P128T2
– P256T1: ピュア MPI実行
• 対象当たりの実行回数：100
72


# Page. 73

![Page Image](https://bcdn.docswell.com/page/KJ42RDK1E1.jpg)

NUMA アフィニティ
• Sparc64 XIfx (FX100) 、および、ARM A64FX (FX1000,
Flow) は NUMA
• (FX100) 2ソケット相当: 16 コア ×2ソケット相当
• (FX1000) 4ソケット相当: 12コア×4ソケット相当
• NUMA アフィニティの詳細 (FX100のみ)
– メモリ割り当て
• “Local allocation”
• plm_ple_memory_allocation_policy=localalloc
– CPU割り当て
• P8 と P16：
plm_ple_numanode_assign_policy=simplex
• P32以上：
plm_ple_numanode_assign_policy=share_band
2023年度 計算科学技術特論A
73


# Page. 74

![Page Image](https://bcdn.docswell.com/page/LE1MQZX57G.jpg)

[Seconds]
200
180
全体時間 （2000 time steps）:FX100
FX100
NX*NY*NZ =
512 x512 x 512 / 8ノード
Comm.
Time
160
Hybrid MPI/OpenMP
FX100: 8ノード (256スレッド)
140
2.55x
120
100
Others
IO
passing_stre
ss
passing_velo
city
Full AT
80
update_vel_
sponge
60
update_vel
40
20
Comp.
Time
update_stre
ss_sponge
オリジナルコード
0
P8T32
P16T16 P32T8
P64T4 P128T2 P256T1
update_stre
ss
74


# Page. 75

![Page Image](https://bcdn.docswell.com/page/GEW3M9NWE2.jpg)

全体時間 （2000 time steps）:FX1000 (「不老」)
[Seconds]
200
180
Comm.
Time
160
140
FX1000 (「不老」)
Others
NX*NY*NZ =
512 x512 x 512 / 8ノード
IO
Hybrid MPI/OpenMP
FX100: 8ノード (384スレッド)
passing_stre
ss
NUMAの影響を
より受けやすい
passing_velo
city
120
100
Full AT
4.06x
80
update_vel_
sponge
60
update_vel
40
20
Comp.
Time
update_stre
ss_sponge
オリジナルコード
0
P8T48
P16T24 P32T12 P64T6
P128T3 P192T2 P384T1
update_stre
ss
75


# Page. 76

![Page Image](https://bcdn.docswell.com/page/47ZVQ932E3.jpg)

そのほかの最適化技法
B/F値の低いカーネル作成
76
2023年度
計算科学技術特論A


# Page. 77

![Page Image](https://bcdn.docswell.com/page/YJ656KZ67V.jpg)

オリジナル実装（ベクトル計算機向き）
速度に関する４次の中央差分スキーム (def_stress)
call ppohFDM_pdiffx3_m4_OAT( VX,DXVX, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
call ppohFDM_pdiffy3_p4_OAT( VX,DYVX, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
call ppohFDM_pdiffz3_p4_OAT( VX,DZVX, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
call ppohFDM_pdiffy3_m4_OAT( VY,DYVY, NXP,NYP,NZP,NXP0,NXP1,NYP0,… )
call ppohFDM_pdiffx3_p4_OAT( VY,DXVY, NXP,NYP,NZP,NXP0,NXP1,NYP0,… )
call ppohFDM_pdiffz3_p4_OAT( VY,DZVY, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
call ppohFDM_pdiffx3_p4_OAT( VZ,DXVZ, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
call ppohFDM_pdiffy3_p4_OAT( VZ,DYVZ, NXP,NYP,NZP,NXP0,NXP1,NYP0,… )
call ppohFDM_pdiffz3_m4_OAT( VZ,DZVZ, NXP,NYP,NZP,NXP0,NXP1,NYP0,…)
if( is_fs .or. is_nearfs ) then
call ppohFDM_bc_vel_deriv( KFSZ,NIFS,NJFS,IFSX,IFSY,IFSZ,JFSX,JFSY,JFSZ )
end if
モデル境界処理
call ppohFDM_update_stress(1, NXP, 1, NYP, 1, NZP)
leap-frogスキームによる陽的時間発展 (update_stress)


# Page. 78

![Page Image](https://bcdn.docswell.com/page/GJ5Y6PW5J4.jpg)

オリジナル実装（ベクトル計算機向き）
subroutine OAT_InstallppohFDMupdate_stress(..)
!$omp parallel do private(i,j,k,RL1,RM1,RM2,RLRM2,DXVX1,DYVY1,DZVZ1,…)
do k = NZ00, NZ01
do j = NY00, NY01
do i = NX00, NX01
RL1 = LAM (I,J,K); RM1 = RIG (I,J,K); RM2 = RM1 + RM1; RLRM2 = RL1+RM2
DXVX1 = DXVX(I,J,K); DYVY1 = DYVY(I,J,K); DZVZ1 = DZVZ(I,J,K)
D3V3 = DXVX1 + DYVY1 + DZVZ1
SXX (I,J,K) = SXX (I,J,K) + (RLRM2*(D3V3)-RM2*(DZVZ1+DYVY1) ) * DT
SYY (I,J,K) = SYY (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DZVZ1) ) * DT
SZZ (I,J,K) = SZZ (I,J,K) + (RLRM2*(D3V3)-RM2*(DXVX1+DYVY1) ) * DT
DXVYDYVX1 = DXVY(I,J,K)+DYVX(I,J,K); DXVZDZVX1 =
DXVZ(I,J,K)+DZVX(I,J,K)
DYVZDZVY1 = DYVZ(I,J,K)+DZVY(I,J,K)
SXY (I,J,K) = SXY (I,J,K) + RM1 * DXVYDYVX1 * DT
SXZ (I,J,K) = SXZ (I,J,K) + RM1 * DXVZDZVX1 * DT
SYZ (I,J,K) = SYZ (I,J,K) + RM1 * DYVZDZVY1 * DT
各手続き呼び出しで
end do
leap-frogスキームに
end do
入出力の配列が多い
よる陽的時間発展
end do
retuen
⇒ B/F値の増加:
(update_stress)
78
end
~1.7


# Page. 79

![Page Image](https://bcdn.docswell.com/page/9E2GM6QGER.jpg)

別実装（スカラ計算機向き）

実装１(IF文あり)


以下の処理をループ中に記載:
1.
速度に関する４次の中央差分スキーム
2.
モデル境界処理
3.
leap-frogスキームによる陽的時間発展
実装２(IF文なし。モデル境界処理のIF文ループが別途必要)
実装１からIF文を取り除くため、ループを再構築。
 計算順序を変更する。しかし、丸め誤差を除く演算結果は同一。
 [主ループ]
1.
速度に関する４次の中央差分スキーム
2.
leap-frogスキームによる陽的時間発展
 [モデル境界処理のループ]
1.
速度に関する４次の中央差分スキーム
2.
モデル境界処理
leap-frogスキームによる陽的時間発展
793.



# Page. 80

![Page Image](https://bcdn.docswell.com/page/D7YVX9WG7M.jpg)

実装１ (スカラ計算機向き)
Stress tensor of Sxx, Syy, Szz
!$omp parallel do private
(i,j,k,RL1,RM1,RM2,RLRM2,DXVX, …)
do k_j=1, (NZ01-NZ00+1)*(NY01-NY00+1)
k=(k_j-1)/(NY01-NY00+1)+NZ00
j=mod((k_j-1),(NY01-NY00+1))+NY00
do i = NX00, NX01
RL1 = LAM (I,J,K); RM1 = RIG (I,J,K)
RM2 = RM1 + RM1; RLRM2 = RL1+RM2
! 4th order diff
DXVX0 = (VX(I,J,K) -VX(I-1,J,K))*C40/dx &amp;
- (VX(I+1,J,K)-VX(I-2,J,K))*C41/dx
DYVY0 = (VY(I,J,K) -VY(I,J-1,K))*C40/dy &amp;
- (VY(I,J+1,K)-VY(I,J-2,K))*C41/dy
DZVZ0 = (VZ(I,J,K) -VZ(I,J,K-1))*C40/dz &amp;
- (VZ(I,J,K+1)-VZ(I,J,K-2))*C41/dz
! truncate_diff_vel
! X dir
if (idx==0) then
if (i==1)then
DXVX0 = ( VX(1,J,K) - 0.0_PN )/ DX
end if
if (i==2) then
DXVX0 = ( VX(2,J,K) - VX(1,J,K) )/ DX
end if
end if
if( idx == IP-1 ) then
if (i==NXP)then
DXVX0 = ( VX(NXP,J,K) - VX(NXP-1,J,K) ) / DX
end if
end if
速度に関する
４次の中央差分
スキーム (def_stress)
80
モデル境界処理
! Y dir
leap-frogスキームに
よる陽的時間発展
(update_stress)
if( idy == 0 ) then ! Shallowmost
if (j==1)then
DYVY0 = ( VY(I,1,K)
- 0.0_PN
)/ DY DYVY1 = DYVY0
DXVX1
= DXVX0;
end if
DZVZ1 = DZVZ0; D3V3 = DXVX1 + DYVY1 +
if (j==2)then
DZVZ1
DYVY0 = ( VY(I,2,K)SXX
- VY(I,1,K)
/ DY (I,J,K) &amp;
(I,J,K) =) SXX
end if
+ (RLRM2*(D3V3)-RM2*(DZVZ1+DYVY1) ) * DT
end if
SYY (I,J,K) = SYY (I,J,K) &amp;
if( idy == JP-1 ) then + (RLRM2*(D3V3)-RM2*(DXVX1+DZVZ1) ) * DT
if (j==NYP)then
SZZ (I,J,K) = SZZ (I,J,K) &amp;
DYVY0 = ( VY(I,NYP,K)
- VY(I,NYP-1,K) )/ DY
+ (RLRM2*(D3V3)-RM2*(DXVX1+DYVY1)
) * DT
end if
end do
end if
end do
! Z dir
!$omp end parallel do
if( idz == 0 ) then ! Shallowmost
if (k==1)then
DZVZ0 = ( VZ(I,J,1) - 0.0_PN ) / DZ
end if
if (k==2)then
DZVZ0 = ( VZ(I,J,2) - VZ(I,J,1) ) / DZ
end if
end if
if( idz == KP-1 ) then
if (k==NZP)then
DZVZ0 = ( VZ(I,J,NZP) - VZ(I,J,NZP-1) )/ DZ
end if
end if
B/F値が
0.4 まで減少
2023年度
計算科学技術特論A
IF文が内在
ーコンパイラに
よる最適化が
難しい


# Page. 81

![Page Image](https://bcdn.docswell.com/page/VENM1L98J8.jpg)

実装2 (IF文なし)
Stress tensor of Sxx, Syy, Szz
!$omp parallel do private(i,j,k,RL1,RM1,･･･)
do k_j=1, (NZ01-NZ00+1)*(NY01-NY00+1)
k=(k_j-1)/(NY01-NY00+1)+NZ00
j=mod((k_j-1),(NY01-NY00+1))+NY00
速度に関する
do i = NX00, NX01
４次の中央差分
RL1 = LAM (I,J,K); RM1 = RIG (I,J,K);
RM2 = RM1 + RM1; RLRM2 = RL1+RM2 スキーム (def_stress)
! 4th order diff (DXVX,DYVY,DZVZ)
DXVX0 = (VX(I,J,K) -VX(I-1,J,K))*C40/dx - (VX(I+1,J,K)-VX(I-2,J,K))*C41/dx
DYVY0 = (VY(I,J,K) -VY(I,J-1,K))*C40/dy - (VY(I,J+1,K)-VY(I,J-2,K))*C41/dy
DZVZ0 = (VZ(I,J,K) -VZ(I,J,K-1))*C40/dz - (VZ(I,J,K+1)-VZ(I,J,K-2))*C41/dz
DXVX1 = DXVX0; DYVY1 = DYVY0;
DZVZ1 = DZVZ0;
D3V3 = DXVX1 + DYVY1 + DZVZ1;
SXX (I,J,K) = SXX (I,J,K) + (RLRM2*(D3V3)-RM2* (DZVZ1+DYVY1) ) * DT
SYY (I,J,K) = SYY (I,J,K) + (RLRM2*(D3V3)-RM2* (DXVX1+DZVZ1) ) * DT
SZZ (I,J,K) = SZZ (I,J,K) + (RLRM2*(D3V3)-RM2* (DXVX1+DYVY1) ) * DT
end do
end do
leap-frogスキームに
!$omp end parallel do
よる陽的時間発展
(update_stress)
B/F値とコンパイラ
最適化の観点から
良いコード


# Page. 82

![Page Image](https://bcdn.docswell.com/page/Y79N542X73.jpg)

実装2 (IF文なし)
! bc vel-derive
if (K==KFSZ(I,J)+1) then
DZVX0 = ( VX(I,J,KFSZ(I,J)+2)-VX(I,J,KFSZ(I,J)+1) )/ DZ
DZVY0 = ( VY(I,J,KFSZ(I,J)+2)-VY(I,J,KFSZ(I,J)+1) )/ DZ
else if (K==KFSZ(I,J)-1) then
DZVX0 = ( VX(I,J,KFSZ(I,J) )-VX(I,J,KFSZ(I,J)-1) )/ DZ
DZVY0 = ( VY(I,J,KFSZ(I,J) )-VY(I,J,KFSZ(I,J)-1) )/ DZ
end if
DXVX1 = DXVX0
DYVY1 = DYVY0
DZVZ1 = DZVZ0
D3V3 = DXVX1 + DYVY1 + DZVZ1
DXVYDYVX1 = DXVY0+DYVX0
DXVZDZVX1 = DXVZ0+DZVX0
DYVZDZVY1 = DYVZ0+DZVY0
if (K==KFSZ(I,J)+1)then
KK=2
else
KK=1
end if
SXX (I,J,K) = SSXX (I,J,KK) &amp;
+ (RLRM2*(D3V3)-RM2*(DZVZ1+DYVY1) ) * D
SYY (I,J,K) = SSYY (I,J,KK) &amp;
+ (RLRM2*(D3V3)-RM2*(DXVX1+DZVZ1) ) * DT
SZZ (I,J,K) = SSZZ (I,J,KK) &amp;
+ (RLRM2*(D3V3)-RM2*(DXVX1+DYVY1) ) * DT
SXY (I,J,K) = SSXY (I,J,KK) + RM1 * DXVYDYVX1 * DT
SXZ (I,J,K) = SSXZ (I,J,KK) + RM1 * DXVZDZVX1 * DT
SYZ (I,J,K) = SSYZ (I,J,KK) + RM1 * DYVZDZVY1 * DT
end do
end do
end do
!$omp end parallel do
モデル境界処理のためのループ
! 2nd replace
if( is_fs .or. is_nearfs ) then
!$omp parallel do private(i,j,k,RL1,RM1,・・・)
do i=NX00,NX01
do j=NY00,NY01
do k = KFSZ(i,j)-1, KFSZ(i,j)+1, 2
RL1 = LAM (I,J,K); RM1 = RIG (I,J,K);
RM2 = RM1 + RM1; RLRM2 = RL1+RM2
! 4th order diff
DXVX0 = (VX(I,J,K) -VX(I-1,J,K))*C40/dx &amp;
- (VX(I+1,J,K)-VX(I-2,J,K))*C41/dx
DYVX0 = (VX(I,J+1,K)-VX(I,J,K) )*C40/dy &amp;
- (VX(I,J+2,K)-VX(I,J-1,K))*C41/dy
DXVY0 = (VY(I+1,J,K)-VY(I ,J,K))*C40/dx &amp;
- (VY(I+2,J,K)-VY(I-1,J,K))*C41/dx
DYVY0 = (VY(I,J,K) -VY(I,J-1,K))*C40/dy &amp;
- (VY(I,J+1,K)-VY(I,J-2,K))*C41/dy
DXVZ0 = (VZ(I+1,J,K)-VZ(I ,J,K))*C40/dx &amp;
- (VZ(I+2,J,K)-VZ(I-1,J,K))*C41/dx
DYVZ0 = (VZ(I,J+1,K)-VZ(I,J,K) )*C40/dy &amp;
- (VZ(I,J+2,K)-VZ(I,J-1,K))*C41/dy
DZVZ0 = (VZ(I,J,K) -VZ(I,J,K-1))*C40/dz &amp;
- (VZ(I,J,K+1)-VZ(I,J,K-2))*C41/dz
速度に関する
４次の中央差分
82
スキーム
(def_stress)
モデル
境界
処理
leap-frogスキーム
による
陽的時間発展
(update_stress)


# Page. 83

![Page Image](https://bcdn.docswell.com/page/G78KYQ59JD.jpg)

速度向上（update_stress）
(a)Xeon Phi (8 Nodes)
[速度向上]
4
3
2.75
2 1.56
1
0
83
最大
オリジナル
ループ分割・消滅
実装変更
2.42
1.40
2.18
1.33
4.2x
2.31
1.34
4.21
2.00
1.23
1.92
1.93
1.25


# Page. 84

![Page Image](https://bcdn.docswell.com/page/L7LZ4XYMER.jpg)

そのほかの最適化技法
データ構造の変換：AoS、SoA
84
2023年度
計算科学技術特論A


# Page. 85

![Page Image](https://bcdn.docswell.com/page/4EMVKLN6EW.jpg)

データ構造変換：AoS vs. SoA

構造体配列：Array of Structures（AoS）
配列

構造体
x0 y0 z0
構造体
x1 y1 z1
構造体
x2 y2 z2
…
構造体
xn yn zn
 利点
 自然なデータ構造で扱いやすい
 欠点
 x0→x1→ｘ2→・・・→xnと
アクセスするとキャッシュミスが生じ、
性能劣化する
配列構造体：Structure of Arrays (SoA)
構造体
85
x配列
x0 x1 x2 x3
…
xn
y配列
y0 y1 y2 y3
…
yn
z配列
z0 z1 z2 z3
…
zn
2023年度
計算科学技術特論A


# Page. 86

![Page Image](https://bcdn.docswell.com/page/PER2XKDN79.jpg)

Modylasにおけるデータ構造変換による
高速化の実例
表
配列形状変更によるP2P(code1)演算時間の削減効果. 単位はms.
1ITO,
-xHost -qopt-zmm-usage=high指定. 2Oakforest-PACS, -axMIC-AVX512指定.
Skylake-SP1
従来
(AoS)
2.140
KNL2
8.582
SoA
実行時間
削減率
1.976
8%
6.048
30%
(source )Modylas:
[1] www.modylas.org
[2] Yoshimichi Ando, et.al., “MODYLAS: A Highly Parallelized General-Purpose Molecular Dynamics Simulation Program for
Large-Scale Systems with Long-Range Forces Calculated by Fast Multipole Method (FMM) and Highly Scalable Fine-Grained
New Parallel Processing Algorithms”, J. Chem. Theo. Comp.,9, 3201-3209 (2013)
[3] Yoshimichi Andoh, et.al., “A thread-level parallelization of pairwise additive potential and force calculations suitable
for current many-core architectures”, The Journal of Supercomputing, 74 (6), pp. 2449-2469 (2018)
86
2023年度
計算科学技術特論A


# Page. 87

![Page Image](https://bcdn.docswell.com/page/P7XD6L18JX.jpg)

AoS からSoAのデータ変換の注意点



実際のプログラムでは、元のデータ構造がAoSになって
いることが多く、プログラム全体をSoAに書き換えること
ができないことが多い
そのため、以下の手順を取ることが多い
1. 対象の計算に入る前にSoAの配列を確保
2. AoSのデータをSoAに変換（コピー）
3. SoAで計算
4. SoAのデータをAoSに変換（コピー）
１により２倍のメモリ量が必要、２、４により、コピー時間
が必要。このデメリットに対して、３の計算時間の
高速化で元が取れる場合のみ有効
87
2023年度
計算科学技術特論A


# Page. 88

![Page Image](https://bcdn.docswell.com/page/37KX1L2LJD.jpg)

通信最適化の方法
88
2023年度
計算科学技術特論A


# Page. 89

![Page Image](https://bcdn.docswell.com/page/LJ35Q34675.jpg)

メッセージサイズと通信回数
領域②
メッセージサイズに比例して、実行時間が
増えていく領域
通信時間[秒]
1 / 傾き係数[秒/バイト]
= メモリバンド幅 [バイト/秒]
領域①
メッセージサイズに
依存せず、ほぼ
一定時間の領域
通信
立ち上がり
時間
=
通信
レイテンシ
[秒]
０
89
領域②の通信時間の計算式
通信時間 = 通信レイテンシ２ ＋
傾き係数 × メッセージサイズ
通信レイテンシ２ [秒]
数百バイト
2023年度
メッセージサイズ[バイト]
計算科学技術特論A


# Page. 90

![Page Image](https://bcdn.docswell.com/page/8JDD94QMJG.jpg)

通信最適化時に注意すること（その１）

自分のアプリケーションの通信パターンについて、
以下の観点を知らないと通信の最適化ができない



領域①の場合



＜領域①＞＜領域②＞のどちらになるのか
通信の頻度（回数）はどれほどか
「通信レイテンシ」が実行時間のほとんど
通信回数を削減する
 細切れに送っているデータをまとめて１回にする、など
領域②の場合



90
「メッセージ転送時間」が実行時間のほとんど
メッセージサイズを削減する
冗長計算をして計算量を増やしてでもメッセージサイズを削減する、など
2023年度
計算科学技術特論A


# Page. 91

![Page Image](https://bcdn.docswell.com/page/VEPXZMLQJ8.jpg)

領域①となる通信の例



内積演算のためのリダクション(MPI_Allreduce)などの送信データ
は倍精度１個分（8バイト）
8バイトの規模だと、数個分を同時にMPI_Allreduceする時間と、
１個分をMPI_Alｌreduceをする時間は、ほぼ同じ時間となる
 ⇒複数回分の内積演算を一度に行うと高速化される可能性あり
例）連立一次方程式の反復解法CG法中の内積演算
 通常の実装だと、１反復に３回の内積演算がある
 このため、内積部分は通信レイテンシ律速となる
 k反復を１度に行えば、内積に関する通信回数は１/k回に削減
 ただし、単純な方法では、丸め誤差の影響で収束しない。
 通信回避CG法（Communication Avoiding CG, CACG）として
現在活発に研究されている。
91
2023年度
計算科学技術特論A


# Page. 92

![Page Image](https://bcdn.docswell.com/page/27VLW9QPEQ.jpg)

通信最適化時に注意すること（その２）

「同期点」を減らすことも高速化につながる



MPI関数の「ノン・ブロッキング関数」を使う
例： ブロッキング関数 MPI_SEND()
→ ノン・ブロッキング関数 MPI_ISEND()
通信と演算を同時に行うようにする。
ランク0
ランク1
計算
計算
send
send
計算
受信待
recv
recv
計算
send
計算
受信待
同期点
受信待
recv
計算
…
…
ノン・ブロッキング関数の利用
ランク0
ランク1
92
計算
計算
isend
irecv
isend
計算
計算
2023年度
計算
irecv
isend
計算
計算科学技術特論A
…
irecv …
高速化


# Page. 93

![Page Image](https://bcdn.docswell.com/page/5JGXDZWQ7L.jpg)

非同期通信：
Isend、Irecv、永続的通信関数
93
2023年度
計算科学技術特論A


# Page. 94

![Page Image](https://bcdn.docswell.com/page/47Q81L3WEP.jpg)

ブロッキング通信で効率の悪い例

プロセス0が必要なデータを持っている場合
連続するsendで、効率の悪い受信待ち時間が多発
プロセス0
プロセス1
計算
プロセス２
計算
プロセス３
計算
…
send
計算
受信待
send
recv
受信待
send
受信待
計算
recv
…
計算
次の反復での同期待ち
計算
recv
次の反復での同期待ち
計算
次の反復での
同期待
同期待ち
…
次の
反復での
同期点
94
2023年度
計算科学技術特論A


# Page. 95

![Page Image](https://bcdn.docswell.com/page/KE42RD1171.jpg)

１対１通信に対するMPI用語
ブロッキング？ノンブロッキング？
95
2023年度
計算科学技術特論A


# Page. 96

![Page Image](https://bcdn.docswell.com/page/L71MQZG5JG.jpg)

ブロッキング、ノンブロッキング
1. ブロッキング

送信／受信側のバッファ領域にメッセージが
格納され、受信／送信側のバッファ領域が
自由にアクセス・上書きできるまで、
呼び出しが戻らない

バッファ領域上のデータの一貫性を保障
ノンブロッキング
2.
96

送信／受信側のバッファ領域のデータを保障せず
すぐに呼び出しが戻る

バッファ領域上のデータの一貫性を保障せず
 一貫性の保証はユーザの責任
2023年度
計算科学技術特論A


# Page. 97

![Page Image](https://bcdn.docswell.com/page/G7W3M9KW72.jpg)

ローカル、ノンローカル

ローカル



ノンローカル


97
手続きの完了が、それを実行しているプロセス
のみに依存する。
ほかのユーザプロセスとの通信を必要としない
処理。
操作を完了するために、別のプロセスでの何らか
のMPI手続きの実行が必要かもしれない。
別のユーザプロセスとの通信を必要とするかもし
れない処理。
2023年度
計算科学技術特論A


# Page. 98

![Page Image](https://bcdn.docswell.com/page/4JZVQ9Z273.jpg)

通信モード（送信発行時の場合）
標準通信モード （ノンローカル） ：デフォルト
1.
送出メッセージのバッファリングはMPIに任せる。



バッファ通信モード （ローカル）
2.
必ずバッファリングする。バッファ領域がないときはエラー。

同期通信モード （ノンローカル）
3.
バッファ領域が再利用でき、かつ、対応する受信／送信が
開始されるまで待つ。

レディ通信モード
4.
（処理自体はローカル）
対応する受信が既に発行されている場合のみ実行できる。
それ以外はエラー。


98
バッファリングされるとき：相手の受信起動前に送信を完了可能；
バッファリングされないとき：送信が完全終了するまで待機；
ハンドシェーク処理を無くせるため、高い性能を発揮する。
2023年度
計算科学技術特論A


# Page. 99

![Page Image](https://bcdn.docswell.com/page/YE656K86JV.jpg)

実例－MPI_Send
 MPI_Send関数


99
ブロッキング
標準通信モード（ノンローカル）
 バッファ領域が安全な状態になるまで戻らない
 バッファ領域がとれる場合：
メッセージがバッファリングされる。対応する受信が
起動する前に、送信を完了できる。
 バッファ領域がとれない場合：
対応する受信が発行されて、かつ、メッセージが
受信側に完全にコピーされるまで、送信処理を
完了できない。
2023年度
計算科学技術特論A


# Page. 100

![Page Image](https://bcdn.docswell.com/page/GE5Y6P95E4.jpg)

非同期通信関数

ierr = MPI_Isend(sendbuf, icount, datatype,
idest, itag, icomm, irequest);





sendbuf : 送信領域の先頭番地を指定する
icount : 整数型。送信領域のデータ要素数を指定する
datatype : 整数型。送信領域のデータの型を指定する
idest : 整数型。送信したいPEのicomm 内でのランクを
指定する
itag : 整数型。受信したいメッセージに付けられたタグ
の値を指定する
100
2023年度
計算科学技術特論A


# Page. 101

![Page Image](https://bcdn.docswell.com/page/972GM62G7R.jpg)

非同期通信関数
 icomm
: 整数型。PE集団を認識する番号
であるコミュニケータを指定する。
通常ではMPI_COMM_WORLD を指定
すればよい。
 irequest : MPI_Request型（整数型の配列）。
送信を要求したメッセージにつけられた
識別子が戻る。
 ierr : 整数型。エラーコードが入る。
101
2023年度
計算科学技術特論A


# Page. 102

![Page Image](https://bcdn.docswell.com/page/DJYVX9YGJM.jpg)

同期待ち関数

ierr = MPI_Wait(irequest, istatus);


irequest : MPI_Request型（整数型配列）。
送信を要求したメッセージにつけられた識別子。
istatus : MPI_Status型（整数型配列）。
受信状況に関する情報が入る。


102
要素数がMPI_STATUS_SIZEの整数配列を宣言して
指定する。
受信したメッセージの送信元のランクが
istatus[MPI_SOURCE] 、タグがistatus[MPI_TAG] に
代入される。
2023年度
計算科学技術特論A


# Page. 103

![Page Image](https://bcdn.docswell.com/page/V7NM1LP8E8.jpg)

実例－MPI_Isend
 MPI_Isend関数


ノンブロッキング
標準通信モード（ノンローカル）
 通信バッファ領域の状態にかかわらず戻る
 バッファ領域がとれる場合は、メッセージが
バッファリングされ、対応する受信が起動する前に、
送信処理が完了できる
 バッファ領域がとれない場合は、対応する受信が
発行され、メッセージが受信側に完全にコピーされる
まで、送信処理が完了できない

103
MPI_Wait関数が呼ばれた場合の振舞いと理解すべき。
2023年度
計算科学技術特論A


# Page. 104

![Page Image](https://bcdn.docswell.com/page/YJ9N543XJ3.jpg)

注意点
 以下のように解釈してください：
 MPI_Send関数
 関数中にMPI_Wait関数が入っている；
 MPI_Isend関数
 関数中にMPI_Wait関数が入っていない；
 かつ、すぐにユーザプログラム戻る；
104
2023年度
計算科学技術特論A


# Page. 105

![Page Image](https://bcdn.docswell.com/page/GJ8KYQM9ED.jpg)

並列化の注意（MPI_Send、MPI_Recv）

全員がMPI_Sendを先に発行すると、その場所で処理が
止まる。(cf. 標準通信モードを考慮)
（正確には、動いたり、動かなかったり、する）
 MPI_Sendの処理中で、場合により、バッファ領域がなくなる。
 バッファ領域が空くまで待つ（スピンウェイトする）。
 しかし、送信側バッファ領域不足から、永遠に空かない。

これを回避するためには、例えば以下の実装を行う。

ランク番号が２で割り切れるプロセス：



それ以外：


105
MPI_Send();
MPI_Recv();
MPI_Recv();
MPI_Send();
それぞれに対応
2023年度
計算科学技術特論A


# Page. 106

![Page Image](https://bcdn.docswell.com/page/LJLZ4X3M7R.jpg)

非同期通信 TIPS

メッセージを完全に受け取ることなく、
受信したメッセージの種類を確認したい

送信メッセージの種類により、受信方式を
変えたい場合
MPI_Probe 関数 （ブロッキング）
 MPI_Iprobe 関数 （ノンブロッキング）
 MPI_Cancel 関数 （ノンブロッキング、
ローカル）

106
2023年度
計算科学技術特論A


# Page. 107

![Page Image](https://bcdn.docswell.com/page/47MVKLD67W.jpg)

MPI_Probe 関数

ierr = MPI_Probe(isource, itag, icomm,
istatus) ;





107
isource: 整数型。送信元のランク。
 MPI_ANY_SOURCE (整数型)も指定可能
itag: 整数型。タグ値。
 MPI_ANY_TAG (整数型) も指定可能
icomm: 整数型。コミュニケータ。
istatus： ステータスオブジェクト。
isource, itagに指定されたものがある場合のみ戻る
2023年度
計算科学技術特論A


# Page. 108

![Page Image](https://bcdn.docswell.com/page/P7R2XK4NJ9.jpg)

MPI_Iprobe関数

ierr = MPI_Iprobe(isource, itag, icomm,
iflag, istatus) ;





108
isource: 整数型。送信元のランク。
 MPI_ANY_SOURCE (整数型) も指定可能。
itag: 整数型。タグ値。
 MPI_ANY_TAG (整数型) も指定可能。
icomm: 整数型。コミュニケータ。
iflag: 論理型。 isource, itagに指定されたものが
あった場合はtrueを返す。
istatus： ステータスオブジェクト。
2023年度
計算科学技術特論A


# Page. 109

![Page Image](https://bcdn.docswell.com/page/PJXD6L28EX.jpg)

MPI_Cancel 関数
 ierr
= MPI_Cancel(irequest);
 irequest:


整数型。通信要求(ハンドル)
目的とする通信が実際に取り消される以前に、
可能な限りすばやく戻る。
取消しを選択するため、MPI_Request_free関数、
MPI_Wait関数、又は MPI_Test関数
（または任意の対応する操作）の呼出しを利用して
完了されている必要がある。
109
2023年度
計算科学技術特論A


# Page. 110

![Page Image](https://bcdn.docswell.com/page/3JKX1LMLED.jpg)

ノン・ブロッキング通信例（C言語）
if (myid == 0) {
…
for (i=１; i&lt;numprocs; i++) {
ierr = MPI_Isend( &amp;a[0], N, MPI_DOUBLE, i,
i_loop, MPI_COMM_WORLD, &amp;irequest[i] );
ランク０のプロセスは、
ランク１~numprocs-1までのプロセス
に対して、ノンブロッキング通信を
用いて、長さNのDouble型配列
データを送信
}
} else {
ierr = MPI_Recv( &amp;a[0], N, MPI_DOUBLE, 0, i_loop,
MPI_COMM_WORLD, &amp;istatus );
}
プロセス0は、recvを
a[ ]を使った計算処理;
待たず計算を開始
if (myid == 0) {
for (i=１; i&lt;numprocs; i++) {
ierr = MPI_Wait(&amp;irequest[i], &amp;istatus);
}
ランク１~numprocs-1までの
プロセスは、ランク０からの
受信待ち。
ランク０のPEは、
ランク１~numprocs-1までのプロセス
に対するそれぞれの送信に対し、
それぞれが受信完了するまで
ビジーウェイト（スピンウェイト）
する。
}
110
2023年度
計算科学技術特論A


# Page. 111

![Page Image](https://bcdn.docswell.com/page/LE35Q3Y6J5.jpg)

ノン・ブロッキング通信の例
（Fortran言語）
if (myid .eq. 0) then
…
do i=１, numprocs - １
call MPI_ISEND( a, N, MPI_DOUBLE_PRECISION,
i, i_loop, MPI_COMM_WORLD, irequest, ierr )
enddo
else
call MPI_RECV( a, N, MPI_DOUBLE,_PRECISION ,
0, i_loop, MPI_COMM_WORLD, istatus, ierr )
endif
a( )を使った計算処理
if (myid .eq. 0) then
プロセス0は、recvを
待たず計算を開始
do i=１, numprocs - １
call MPI_WAIT(irequest(i), istatus, ierr )
enddo
ランク０のプロセスは、
ランク１~numprocs-1までの
プロセスに対して、ノンブロッキング
通信を用いて、長さNの
DOUBLE PRECISION型配列
データを送信
ランク１~numprocs-1までの
プロセスは、
ランク０からの受信待ち。
ランク０のプロセスは、
ランク１~numprocs-1までの
プロセスに対するそれぞれの送信
に対し、それぞれが受信完了
するまでビジーウェイト
（スピンウェイト）する。
endif
111
2023年度
計算科学技術特論A


# Page. 112

![Page Image](https://bcdn.docswell.com/page/8EDD945MEG.jpg)

ノン・ブロッキング通信による改善

プロセス0が必要なデータを持っている場合
連続するsendにおける受信待ち時間を
ノン・ブロッキング通信で削減
プロセス0
プロセス1
計算
プロセス２
計算
プロセス３
計算
…
112
send send send
計算
recv
…
計算
受信待ちを、MPI_Waitで
計算の後に行うように変更
計算
受信待
次の反復で
同期待ち
の同期待ち
次の反復での同期待ち
recv
計算
recv
…
次の反復での同期待ち
計算
次の
反復での
同期点
2023年度
計算科学技術特論A


# Page. 113

![Page Image](https://bcdn.docswell.com/page/V7PXZMGQE8.jpg)

永続的通信（その１）


ノン・ブロッキング通信は、MPI_ISENDの実装が、
MPI_ISENDを呼ばれた時点で本当に通信を開始する
実装になっていないと意味がない。
ところが、MPIの実装によっては、MPI_WAITが呼ばれる
まで、MPI_ISENDの通信を開始しない実装がされている
ことがある。


この場合には、ノン・ブロッキング通信の効果が全くない。
永続的通信（Persistent Communication）を利用すると、
MPIライブラリの実装に依存し、ノン・ブロッキング通信の
効果が期待できる場合がある。

永続的通信は、MPI-１からの仕様（たいていのMPIで使える）

113
しかし、通信と演算がオーバラップできる実装になっているかは別問題
2023年度
計算科学技術特論A


# Page. 114

![Page Image](https://bcdn.docswell.com/page/2JVLW96P7Q.jpg)

永続的通信（その２）

永続的通信の利用法
1.
2.
3.

MPI_SEND_INIT関数で通信情報を設定しておくと、
MPI_START時に通信情報の設定が行われない


通信を利用するループ等に入る前に１度、通信相手先を
設定する初期化関数を呼ぶ
その後、SENDをする箇所にMPI_START関数を書く
真の同期ポイントに使う関数(MPI_WAIT等)は、ISENDと同じ
ものを使う
同じ通信相手に何度でもデータを送る場合、通常の
ノン・ブロッキング通信に対し、同等以上の性能が出ると期待
適用例


114
領域分割に基づく陽解法
陰解法のうち反復解法を使っている数値解法
2023年度
計算科学技術特論A


# Page. 115

![Page Image](https://bcdn.docswell.com/page/5EGXDZNQJL.jpg)

永続的通信の実装例（C言語）
MPI_Status istatus;
メインループに入る前に、
MPI_Request irequest;
送信データの相手先情報を
…
初期化する
if (myid == 0) {
for (i=1; i&lt;numprocs; i++) {
ierr = MPI_Send_init (a, N, MPI_DOUBLE_PRECISION, i,
0, MPI_COMM_WORLD, irequest );
}
}
…
ここで、データを送る
if (myid == 0) {
for (i=1; i&lt;numprocs; i++) {
ierr = MPI_Start ( irequest );
}
}
/* 以降は、Isendの例と同じ */
115
2023年度
計算科学技術特論A


# Page. 116

![Page Image](https://bcdn.docswell.com/page/4JQ81LQW7P.jpg)

永続的通信の実装例（Fortran言語）
integer istatus(MPI_STATUS_SIZE)
メインループに入る前に、
integer irequest(0:MAX_RANK_SIZE)
送信データの相手先情報を
…
初期化する
if (myid .eq. 0) then
do i=1, numprocs-1
call MPI_SEND_INIT (a, N, MPI_DOUBLE_PRECISION, i,
0, MPI_COMM_WORLD, irequest(i), ierr)
enddo
endif
…
ここで、データを送る
if (myid .eq. 0) then
do i=1, numprocs-1
call MPI_START (irequest, ierr )
enddo
endif
/* 以降は、ISENDの例と同じ */
116
2023年度
計算科学技術特論A


# Page. 117

![Page Image](https://bcdn.docswell.com/page/K742RDN1J1.jpg)

大規模学習への展開
117
2023年度
計算科学技術特論A


# Page. 118

![Page Image](https://bcdn.docswell.com/page/LJ1MQZ55EG.jpg)

大規模GPUクラスタによる
機械学習ハイパーパラメタ最適化例



機械学習時に事前設定する
パラメタ（ハイパーパラメタ）の最適化
パラメタを変動させて、多数の機械学習を実行する
必要 →膨大な実行時間
機械学習の速さではなく（当然、実行時間も
速いほうが良いが）、CNNの精度向上が目的


一度最適なパラメタでCNNを生成したら、それを
使いまわすような利用法
パラメタチューニングとみなせる
藤家、田部田、藤井、田中（工学院大）、加藤 （東女大）、大島、片桐（名大）「GPUクラスタ
を用いて並列化した自動チューニングの機械学習プログラムへの適用と安定性の検証 」
研究報告ハイパフォーマンスコンピューティング（HPC）、vol. 2021-HPC-178、No. 16、
pp.1 – 8 (2021-03-08)
118
2023年度
計算科学技術特論A


# Page. 119

![Page Image](https://bcdn.docswell.com/page/GJW3M95WJ2.jpg)

工学院大学 田中研究室提供
自動チューニング手法
• 性能パラメータ空間
► 複数の性能パラメータが取り得る値
からなる離散空間
• 多次元の性能パラメータ空間において
方向探索と一次元探索を繰り返す
► 方向探索
•
現在の最良点の各方向の隣接点を調べ，
最も性能の良い点を見つける
開始
初期点を選択
方向探索
一次元探索
探索する軸を追加
► 一次元探索
•
方向探索で見つけた方向にある点を全て調べ，
最も性能の良い点を新たな最良点とする
最良点が
変わらない
No
Yes
• 点の調べ方
► 性能パラメータの組合せを設定し1回実測
全軸
探索終了
No
Yes
終了
望月大義，藤井昭宏，田中輝雄，ソフトウェア自動チューニングにおける複数同時性能
パラメタ探索手法の提案と評価，情報処理学会論文誌，vol.11，No2，pp.1-16（2018）．
119


# Page. 120

![Page Image](https://bcdn.docswell.com/page/4EZVQ9N2J3.jpg)

工学院大学 田中研究室提供
推定フェーズの並列化
• 複数のジョブを同時に実行
できる環境に向けて並列化
• 複数のジョブで同時に実測
►それぞれのジョブが
別の性能パラメータを実測
• 方向探索と一次元探索を
まとめて行う
開始
初期点を選択
並列化
方向探索
一次元探索
探索する軸を追加
最良値が
同点
Yes
全軸
探索終了
Yes
No
No
終了
2023年度 計算科学技術特論A
120


# Page. 121

![Page Image](https://bcdn.docswell.com/page/Y7656K96EV.jpg)

工学院大学 田中研究室提供
並列時のジョブ実行手順
全ての性能パラメータ
を読み込み
ジョブ
ジョブ投入
ジョブごとの
性能パラメータを
設定
チューニング対象
プログラム実⾏
評価値を送信
ジョブ投入確認
全ジョブ完了まで待機
自動チューニング機構
性能パラメータが
取り得る値の範囲を設定
各ジョブの評価値を
まとめて送信
2023年度 計算科学技術特論A
121


# Page. 122

![Page Image](https://bcdn.docswell.com/page/G75Y6PN574.jpg)

工学院大学 田中研究室提供
機械学習への適用
• 自動チューニングを機械学習プログラムに適用
• 機械学習
►訓練データを使って学習を行い，
学習結果から目的となるタスクをこなす
►多くの場合，学習は複数回行われ，
その度に訓練データがランダムに変化
• ハイパーパラメータ
►学習率やドロップアウト率のような機械学習の挙動を制御する
パラメータ
►経験則で調整するのは難しく，膨大な時間が必要
2023年度 計算科学技術特論A
122


# Page. 123

![Page Image](https://bcdn.docswell.com/page/9J2GM65GJR.jpg)

工学院大学 田中研究室提供
対象プログラム
• 歩行者経路予測用アプリケーション
►ロボットが歩行者に追従するように
歩行者の経路と到着地点を予測
►予測した到達地点と実際の到達地点の誤差である
FDE（Final Displacement Error）[m] で評価(小さい程良い)
►Rina Akabane, Yuka Kato, “Pedestrian Trajectory Prediction Using
Pre-trained Machine Learning Model for Human-Following Mobile
Robot”, IEEE International Conference on Big Data Workshop
(IoTDA 2020), pp.3453-3458 (2020).
• つのステップで経路を予測
1.
2.
訓練データ（予め計測しておいたデータ）を用いて
学習を行い，学習結果から予測器を生成
予測器にロボットが計測した移動軌跡を入力し経路を予測
2023年度 計算科学技術特論A
123


# Page. 124

![Page Image](https://bcdn.docswell.com/page/DEYVX9KGEM.jpg)

問題設定
工学院大学 田中研究室提供
チューニング対象のハイパーパラメータ
• 種類のハイパーパラメータが
それぞれ つの取りうる値を持つ
• ハイパーパラメータの組み合わせ：
ହ
通り
ハイパーパラメータ
特徴
Rnn size
Grad clip
Learning rate
Dropout
Lamba
LSTM (Long-Short Term Memory) において短期記憶の役割を果たす hidden state の大きさ
勾配が大きくならないように修正するための閾値
1回の学習でニューラルネットワーク内の重みやバイアスを更新する量の調整値
過学習を抑えるために，学習時に特定のレイヤーの出力を０に落とす割合
過学習を防ぐためのL2正則化における正則化パラメータ
124


# Page. 125

![Page Image](https://bcdn.docswell.com/page/VJNM1LR878.jpg)

工学院大学 田中研究室提供
対象プログラムとの位置付け
自動チューニングにより，
新しい予測器が得られた
段階で，予測器を更新
スーパーコンピュータで
並列実行
予測モデルの生成
Cyber
World
(クラウド)
データ
変換
歩行者の
位置情報
抽出
自動チューニング機構
・反復一次元探索
モデル更新
(ハイパーパラメータ設定)
訓練データ
性能評価値
(FDE)
機械学習により
予測器を更新
・Social LSTM
クラウドーエッジ連携
Physical 位置情報取得
World
センシング
(エッジ)
・LiDAR
[ロボット]
歩行者移動予測
データ変換
(位置座標)
2023年度 計算科学技術特論A
予測器
歩行者移動予測
125


# Page. 126

![Page Image](https://bcdn.docswell.com/page/YE9N54MXE3.jpg)

工学院大学 田中研究室提供
実行環境
• 名古屋大学のスーパーコンピュータ
Ⅱサブシステム
「不老」
機種名
FUJITSU Server PRIMERGY CX2570 M5
Intel Xeon Gold 6230, 2,560 FP64コア,
20コア, 2.10 - 3.90 GHz × 2 ソケット
NVIDIA Tesla V100 (Volta) SXM2, 2,560 FP64コア,
upto 1,530 MHz × 4 ソケット
CPU
GPU
メモリ
理論演算性能
ノード数，
総コア数
メインメモリ(DDR4 2933 MHz) 384 GiB (32 GiB × 6 枚 × 2 ソケット)
デバイスメモリ(HBM2) 32 GiB × 4 ソケット
倍精度 33.888 TFLOPS
(CPU 1.344 TFLOPS × 2 ソケット, GPU 7.8 TFLOPS × 4 ソケット)
221 ノード，
8,840 CPUコア + 2,263,040 FP64 GPUコア
http://www.icts.nagoya-u.ac.jp/ja/sc/overview.html#type2
126


# Page. 127

![Page Image](https://bcdn.docswell.com/page/GE8KYQL97D.jpg)

工学院大学 田中研究室提供
「不老」を用いたジョブ実行
ログインノード
ジョブ
計算ノード
全ての性能パラメータ
を読み込み
ジョブ投入
ジョブごとの
性能パラメータを
設定
チューニング対象
プログラム実⾏
評価値を送信
ジョブ投入確認
全ジョブ完了まで待機
自動チューニング機構
性能パラメータが
取り得る値の範囲を設定
各ジョブの評価値を
まとめて送信
2023年度 計算科学技術特論A
127


# Page. 128

![Page Image](https://bcdn.docswell.com/page/LELZ4XLMJR.jpg)

並列化によるチューニング時間の短縮実験
工学院大学 田中研究室提供
実験 ：並列実測と逐次実測の比較
・探索の初めの基準点となる
性能パラメータ
・経験則に基づいて現在使用
128


# Page. 129

![Page Image](https://bcdn.docswell.com/page/4JMVKLM6JW.jpg)

超解像の機械学習プログラム
工学院大学 田中研究室提供
• 超解像とは，低解像度画像を高解像度画像に変換する
技術
• 対象プログラム:Dense Deep Back-Projection
Networks (D-DBPN) *1)
• 豊田工業大学の浮田宗伯教授のグループが開発
• 低解像度画像の各辺を2倍,4倍, 8倍の高解像度画像に拡大
• 本研究では4倍への拡大
• 画像の画質を性能評価値(Perceptual Index, RMSE)*2) で評価
•
•
•
•
Perceptual Index(PI), RMSEが小さいほど性能が良い
PI は人が見て綺麗かどうかの指標
RMSEは試験画像と予測画像との誤差を計算
MATLAB でPIRM (Perceptual Image Restoration and Manipulation)
プログラムを実行して性能評価値を計算
*1) M. Haris, G. Shakhnarovich and N. Ukita, &quot;Deep Back-Projection Networks for Super-Resolution,“ IEEE/CVF (2018)
*2) Y. Blau, R. Mechrez, R. Timofte, T. Michaeli, and L. Zelnik-Manor, “The 2018 PIRM challenge on perceptual
image super-resolution,” ECCV(2018)
129


# Page. 130

![Page Image](https://bcdn.docswell.com/page/PJR2XKVNE9.jpg)

超解像プログラムの設定
工学院大学 田中研究室提供
• ハイパーパラメータ
• 4種類: 学習率， w2 (Perceptual損失ウェイト)， w3 (Adversarial 損失ウェイト)，w4(Style損失ウェイト)
• w1 (MSE損失ウェイト)を固定: 10ିଵ
• 組み合わせ13,310(ൌ 10 ൈ 11 ൈ 11 ൈ 11), 61,440(ൌ 15 ൈ 16 ൈ 16 ൈ 16)通り
10ି଺ , 5 ൈ 10ି଺ , 10ିହ , 5 ൈ 10ିହ , 10ିସ , 5 ൈ 10ିସ , 10ିଷ , 5 ൈ 10ିଷ , 10ିଶ , 5 ൈ 10ିଶ
10パターン
w2
0, 10ିଶ , 5 ൈ 10ିଶ , 10ିଵ , 5 ൈ 10ିଵ , 1, 5, 10ଵ , 5 ൈ 10ଵ , 10ଶ , 5 ൈ 10ଶ
11パターン
w3
0, 10ିହ , 5 ൈ 10ିହ , 10ିସ , 5 ൈ 10ିସ , 10ିଷ , 5 ൈ 10ିଷ , 10ିଶ , 5 ൈ 10ିଶ , 10ିଵ , 5 ൈ 10ିଵ
11パターン
w4
0, 10ିଶ , 5 ൈ 10ିଶ , 10ିଵ , 5 ൈ 10ିଵ , 1, 5, 10ଵ , 5 ൈ 10ଵ , 10ଶ , 5 ൈ 10ଶ
11パターン
10-6, 2×10-6, 5×10-6, 10-5, 2×10-5, 5×10-5, 10-4, 2×10-4, 5×10-4, 10-3, 2×10-3, 5×10-3, 10-2, 2×10-2, 5×10-2
15パターン
w2
0, 10-2, 2×10-2, 5×10-2, 10-1, 2×10-1, 5×10-1, 1.0, 2.0, 5.0, 10, 20, 50, 100, 200, 500
16パターン
w3
0, 10-5, 2×10-5 , 5×10-5, 10-4, 2×10-4, 5×10-4, 10-3, 2×10-3, 5×10-3, 10-2, 2×10-2, 5×10-2, 10-1, 2×10-1, 5×10-1
16パターン
w4
0, 10-2, 2×10-2, 5×10-2, 10-1, 2×10-1, 5×10-1, 1.0, 2.0, 5.0, 10, 20, 50, 100,200,500
16パターン
学習率
学習率
130


# Page. 131

![Page Image](https://bcdn.docswell.com/page/PEXD6LZ87X.jpg)

工学院大学 田中研究室提供
ジョブの並列実⾏状況
名大「不老」TypeⅡサブシステムの
１ユーザの最大同時ジョブ実行数：50並列
各軸方向の探索
2軸の探索
3軸の探索
4軸の探索
120
45
127
64
149
100
80
60
40
20
45並列
10並列
50並列
0
1
11
21
31
41
51
61
71
81
91
101
111
121
131
141
151
161
171
181
191
201
211
221
231
241
251
261
271
281
291
301
311
321
331
341
351
361
371
381
391
401
411
421
431
441
実行時間[時間]
60
ジョブ番号
総実行回数
総実行時間
445
110h
131


# Page. 132

![Page Image](https://bcdn.docswell.com/page/3EKX1L8L7D.jpg)

結果画像
工学院大学 田中研究室提供
RMSE&lt;11.5
11.5&lt;RMSE&lt;12.5
12.5&lt;RMSE&lt;16
132


# Page. 133

![Page Image](https://bcdn.docswell.com/page/L735Q326E5.jpg)

レポート課題（その１）

問題レベルを以下に設定
問題のレベルに関する記述：
•L00: きわめて簡単な問題。
•L10： ちょっと考えればわかる問題。
•L20： 標準的な問題。
•L30： 数時間程度必要とする問題。
•L40： 数週間程度必要とする問題。複雑な実装を必要とする。
•L50： 数か月程度必要とする問題。未解決問題を含む。
※L４０以上は、論文を出版するに値する問題。

教科書のサンプルプログラムは以下が利用可能

133
付属のサンプルプログラム全てが利用可能
2023年度
計算科学技術特論A


# Page. 134

![Page Image](https://bcdn.docswell.com/page/87DD94ZM7G.jpg)

レポート課題（その２）
1.
2.
3.
4.
[L5] ＭＰＩにおけるブロッキングは、必ずしも同期でないこと
を説明せよ。
[L１0] MPIにおけるブロッキング、ノンブロッキング、および
通信モードによる分類に対応する関数を調べ、一覧表に
まとめよ。
[L１5] 利用できる並列計算機環境で、ノンブロッキング送信
（MPI_Isend関数）がブロッキング送信（MPI_Send関数）に対
して有効となるメッセージの範囲（N=0～適当な上限）につ
いて調べ、結果を考察せよ。
[L20] MPI_Allreduce関数 の＜限定機能＞版を、ブロッキン
グ送信、およびノンブロッキング送信を用いて実装せよ。さ
らに、その性能を比べてみよ。なお、＜限定機能＞は独自
に設定してよい。
134
2023年度
計算科学技術特論A


# Page. 135

![Page Image](https://bcdn.docswell.com/page/VJPXZMDQ78.jpg)

レポート課題（その３）
5.
6.
7.
8.
9.
[L15] MPI_Reduce関数を実現するRecursive Halving
アルゴリズムについて、その性能を調査せよ。この際、
従来手法も調べて、その手法との比較も行うこと。
[L35] Recursive Halvingアルゴリズムを、ブロッキング送信／
受信、および、ノンブロッキング送信／受信を用いて実装せよ。
また、それらの性能を評価せよ。
[L15] 身近の並列計算機環境で、永続的通信関数の性能を
調べよ。
[L10～] 自分が持っているプログラムに対し、ループ分割、
ループ融合、その他のチューニングを試みよ。
[L10～] 自分が持っているMPIプログラムに対し、ノンブロッキン
グ通信（MPI_Isend, MPI_Irecv）を実装し、性能を評価せよ。
また永続的通信が使えるプログラムの場合は実装して評価せよ。
135
2023年度
計算科学技術特論A


