4.1 核融合プラズマ乱流コード GKV の Post-FX100 性能推定
核融合科学研究所 沼波 政倫
名古屋大学情報基盤センター 片桐 孝洋
富士通株式会社

青木 正樹

4.1.1 は じ め に
核 融 合 プ ラ ズ マ 乱 流 コ ー ド で あ る GKV（ ス ラ イ ド p.3） は 、 主 に 磁 場 閉 じ 込 め型の核融合装置における炉心プラズマの乱流輸送現象を扱うように設計・開発された。磁場閉じ込め装置は、強い磁場によってトーラス状の高温プラズマを閉じ込める。この強い磁場に起因して測度空間を 2 次元に近似した 5 次元位相
空 間 上 （ 3 次 元 実 空 間 +2 次 元 速 度 空 間 ） の プ ラ ズ マ 分 布 関 数 の 発 展 を 記 述 す る
ボルツマン方程式（ジャイロ運動論方程式）が、この磁場閉じ込めプラズマの第 一 原 理 モ デ ル と な る 。GKV コ ー ド で は 、ジ ャ イ ロ 運 動 論 方 程 式 を 5 次 元 格 子 上
で離散化して解き進める。特に、乱流による揺らぎの磁力線垂直方向への波長が平衡分布のスケールよりも十分に短いことを仮定し、一本の磁力線に沿った
局 所 フ ラ ッ ク ス・チ ュ ー ブ 上 で シ ミ ュ レ ー シ ョ ン を 行 う 。5 次 元 位 相 空 間 は 、空間離散化として磁力線垂直方向の 2 次元空間に対してはフーリエ・スペクトル
法 、磁 力 線 平 行 方 向 の 実 空 間 、お よ び 2 次 元 の 速 度 空 間 方 向 に つ い て は 、4 次 ま
た は 5 次 精 度 差 分 法 を 用 い て お り 、 時 間 積 分 に は 4 次 精 度 陽 的 Runge-Kutta-Gi
ll 法 を 採 用 し て い る 。 2 次 元 の 速 度 空 間 に 加 え て 、 電 磁 揺 動 や 極 微 細 電 子 ス ケール乱流を扱うことができるようにするため、磁力線垂直方向の 2 次元実空間
は 非 常 に 高 い 解 像 度 を 必 要 と し 、 大 型 ヘ リ カ ル 装 置 （ LHD） な ど の 複 雑 な 磁 場 形状を有する実験装置での解析も対象とするために、磁力線平行方向にも高解像度を要する。さらには、プラズマを構成する粒子の種類の分だけ、分布関数の
種 類 が 増 加 す る た め 、 最 終 的 に 総 格 子 点 数 と し て 100 億 点 以 上 を 必 要 と す る 極
めて大きな自由度の流体計算となる。「京」による全系計算でようやく、電子と水素イオンのスケールの乱流を同時に扱うことが初めてできるようになったが、さらに質量や電荷の大きなイオンを含んだ現実的なプラズマを扱うためには、計算コストがさらに数十倍〜数百倍に増大するため、ポスト京のような次世代エクサスケール計算機が必要である。京以降で現行の計算機、例えば核融合 科 学 研 究 所 の プ ラ ズ マ シ ミ ュ レ ー タ（ FX100：SPARC64XIfx）か ら 、ポ ス ト 京 、
あ る い は 、 Post-FX100（ A64fx） へ の 移 行 を 考 え た 場 合 、 大 き な 技 術 的 ギ ャ ッ プ
が 存 在 す る 。例 え ば 、SPARC か ら ARM へ の ア ー キ テ ク チ ャ の 変 更 や 、レ ジ ス タ 数が 128 か ら 32 へ 削 減 、 SIMD 幅 の 256bit か ら 512bit へ の 増 加 、 プ ロ セ ッ サ （ CMG） あ た り 演 算 コ ア 数 の 32（ 16） か ら 48（ 12） へ の 変 更 等 で あ る 。 こ の う ち 、
本性能推定ではユーザプログラムの実装方法への影響が大きいと思われる、レジ ス タ 数 の 削 減 に 関 し て 、FX100 上 で 用 意 さ れ た ポ ス ト FX100 向 け チ ュ ー ニ ン グ

用 コ ン パ イ ラ と Post-FX100 性 能 推 定 ツ ー ル を 用 い て 、Post-FX100 に 対 す る GKV
コードの性能推定を実施した。

4.1.2 プ ロ グ ラ ム 概 要
GKV は 、プ ラ ズ マ 分 布 関 数 の 熱 平 衡 部 分 と 摂 動 部 分 に 分 け た 場 合 の 摂 動 部 分 の時 間 発 展 を 追 う Delta-f モ デ ル に 基 づ き 、磁 場 閉 じ 込 め 炉 心 プ ラ ズ マ に 対 し て 、一本の磁力線に沿った局所領域（フラックス・チューブ）注目し、この局所領域を計算対象としてプラズマの乱流およびプラズマ分布関数の時間発展を追跡する。大域的な手法に基づく乱流コードと違い、平衡分布の発展は記述することができないが、電磁揺動や極微細電子スケール乱流を定量的に扱うことができ る た め 、 こ れ ま で に 、 大 型 ヘ リ カ ル 装 置 （ LHD） や JT-60U な ど 、 実 際 の 大 型
プラズマ実験装置での放電プラズマを対象にした実験結果の再現や予測解析が行われてきた。前述の通り、空間離散化としては、磁力線垂直方向の 2 次元空間 に 対 し て フ ー リ エ・ス ペ ク ト ル 法 、磁 力 線 平 行 方 向 の 1 次 元 実 空 間 、お よ び 2
次 元 の 速 度 空 間 方 向 に つ い て は 、4 次 ま た は 5 次 精 度 差 分 法 を 採 用 し 、時 間 積 分に は 4 次 精 度 陽 的 Runge-Kutta-Gill 法 を 採 用 し て い る 。 GKV で は 、 上 述 の 5 次
元 位 相 空 間 に 加 え て 、プ ラ ズ マ 粒 子 種 を 区 別 す る た め に 6 次 元 配 列 を 用 意 し 、M
PI を 利 用 し た 5 次 元 領 域 分 割 を 行 う 。 こ の 際 、 並 列 2 次 元 FFT の た め の 転 置 通信、差分演算のための 1 対 1 通信、速度空間・粒子種積分のための総和通信とい っ た MPI 通 信 を 利 用 し て い る 。 こ れ ま で に 、 「 京 」 に 対 す る 最 適 化 と し て 、 T
ofu ネ ッ ト ワ ー ク に お け る プ ロ セ ス 配 置 最 適 化 、 OpenMP を 利 用 し た 通 信 演 算 オ
ー バ ー ラ ッ プ の 実 装 等 の 最 適 化 処 理 が 施 さ れ て お り 、種 々 の 効 率 化 の 結 果 、「京」
全 系 規 模 に 渡 る 約 60 万 コ ア の 強 ス ケ ー リ ン グ （ 演 算 性 能 786.4TFlops、 理 論 ピ
ー ク 比 8.29%、 実 行 並 列 化 率 99.99994%） を 達 成 し 、 前 述 の よ う に 「 京 」 の 全 系
計 算 で 電 子 ス ケ ー ル 乱 流 と イ オ ン ス ケ ー ル 乱 流 の 同 時 計 算 が 実 現 し た 。ま た 、F
X100 に お け る ア シ ス タ ン ト ・ コ ア を 用 い て 、 通 信 と 演 算 の オ ー バ ー ラ ッ プ 最 適
化 も 行 っ て い る 。 1CMG あ た り 16 演 算 コ ア + 1 ア シ ス タ ン ト ・ コ ア を 有 す る FX
100 の 特 性 を 反 映 し て 、ア シ ス タ ン ト・コ ア で 通 信 処 理 を 行 っ た 場 合 、演 算 コ アの マ ス タ ー ス レ ッ ド が 通 信 に 参 加 す る 場 合 と 比 較 し て 、 6.25%（ =1/16） 近 く の
高 速 化 を 達 成 し て い る （ ス ラ イ ド p.4〜 5 参 照 ） 。

4.1.3

ポ ス ト FX100 向 け チ ュ ー ニ ン グ 用 コ ン パ イ ラ に よ る コ ス ト 分 析

GKV コ ー ド の 性 能 推 定 に 向 け て 、 GKV カ ー ネ ル コ ー ド を 対 象 に 、 現 行 の FX100の １ ノ ー ド を 用 い た 単 体 コ ス ト 分 析 を 行 っ た 。 コ ン パ イ ラ は ポ ス ト FX100 向 け
チ ュ ー ニ ン グ 用 コ ン パ イ ラ を 用 い 、 レ ジ ス タ 数 128（ 以 降 、 128reg と 記 載 ） での コ ン パ イ ル で の オ プ シ ョ ン は 、 「 -Kfast -Kparallel -Kocl -Cpp -X9」 で あり 、 レ ジ ス タ 数 32（ 以 降 、 32reg と 記 載 ） で の コ ン パ イ ル ・ オ プ シ ョ ン は 128r
eg と 同 様 に ソ フ ト ウ ェ ア ・ パ イ プ ラ イ ン を 動 作 さ せ る た め 、 上 に 加 え て 「 -Kno

unroll」を 追 加 指 定 し た 。ス ラ イ ド p.7 に 、測 定 条 件 と コ ス ト 分 布 の 例 を 示 す 。
こ こ で は 、１ ノ ー ド １ プ ロ セ ス（ MPI 並 列 な し ）を 前 提 と し て 、取 り 扱 う 問 題 サ
イ ズ を 、 通 常 の １ MPI プ ロ セ ス あ た り に 相 当 す る （ nx,ny,nz,nv,nm） =（ 85,64,
24,8,48） と し た （ 粒 子 種 数 は 1） 。 128reg で の コ ン パ イ ラ を 用 い て 測 定 し た G
KV コ ー ド の コ ス ト 分 布 を 、ス ラ イ ド p.7 下 部 に 示 す 。こ の コ ス ト 分 析 か ら 、「 literm_k」 お よ び 「 exb_realspcal」 の サ ブ ル ー チ ン が そ れ ぞ れ 全 体 の 約 53.7%、
約 13.6%を 占 め て い る た め 、以 降 の コ ス ト 評 価 は 、こ の 2 つ に お け る ル ー プ を 対象にする。
前 者 の サ ブ ル ー チ ン 「 literm_k」 で は 、 4 次 精 度 陽 的 Runge-Kutta-Gill
法による時間積分を行う際の方程式の右辺における線形項計算の一部に相当す
る も の で あ り 、 （ nx,ny,nz,nv,nm） に 関 す る 5 次 元 ル ー プ か ら 構 成 さ れ て い る（ ス ラ イ ド p.9 参 照 ） 。 128reg コ ン パ イ ラ と 32reg コ ン パ イ ラ で の 実 行 時 間 の比 較（ ス ラ イ ド p.10）で は 、128reg か ら 32reg へ の コ ン パ イ ラ の 変 更 に よ る 性
能 の 劣 化 は ほ と ん ど 見 ら れ て い な い 。 ま た 、 詳 細 PA 情 報 を 用 い た 比 較 （ ス ラ イ
ド p.11） で は 、 128reg コ ン パ イ ラ 、 32reg コ ン パ イ ラ と も に 、 メ モ リ ビ ジ ー 率が 高 く （ 96%） 、 メ モ リ バ ン ド 幅 ネ ッ ク の ル ー プ で あ る こ と が 分 か る 。 一 方 、 L1
D キ ャ ッ シ ュ の ミ ス 率 は 、 理 論 値 （ 3.125%） 付 近 で あ っ た 。
後 者 の サ ブ ル ー チ ン 「 exb_realspcal」 で は 、 時 間 積 分 を 行 う 際 の 方 程 式 の 右
辺 に お け る 非 線 形 項 計 算 の 一 部 に 相 当 し て お り 、（ nx,ny,nv,nm） に 関 す る 4 次
元 ル ー プ 計 算 か ら 構 成 さ れ て い る（ ス ラ イ ド p.13 参 照 ）。128reg コ ン パ イ ラ での 実 行 時 間 と 比 較 し て 、32reg コ ン パ イ ラ の 場 合 は 、4%程 度 の 短 縮 と な っ て お り（ ス ラ イ ド p.14） 、 「 literm_k」 の 場 合 と 同 様 に 、 32reg コ ン パ イ ラ に よ る 性
能 劣 化 は ほ ぼ 見 ら れ て い な い 。詳 細 PA 情 報 に よ る 比 較（ ス ラ イ ド p.15）に お いて も 、 128reg お よ び 32reg コ ン パ イ ラ と も に メ モ リ ビ ジ ー 率 が 高 く 、 メ モ リ バ
ン ド 幅 ネ ッ ク の ル ー プ で あ る こ と が 分 か る 。一 方 、L1D ミ ス 率 が そ れ ぞ れ 5.12%、
5.24%と な っ て お り 、 理 論 値 （ 3.15%） よ り も 高 い た め 、 ス ラ ッ シ ン グ が 発 生 して い る と 考 え ら れ る 。こ の 対 応 と し て 、5 節 で 述 べ る ル ー プ 分 割 に よ る チ ュ ー ニ
ングを試みた。

4.1.4

Post-FX100 に お け る 性 能 の 推 定

GKV カ ー ネ ル コ ー ド の Post-FX100 に お け る 性 能 を 推 定 す る に あ た り 、 前 節 のコスト分析で対象にした 2 つのサブルーチンについて詳細な評価を行う。性能
予 測 で は 、 ス ラ イ ド p.17 に 示 す よ う に 、 FX100 に お い て 32reg コ ン パ イ ラ を 利用し、詳細プロファイラの実行でＰＡ情報を取得する。そこから、性能予測ツールを用いて、ハード、命令セットおよびコンパイラによる性能向上効果を算出し、性能予測値を出力する。
も っ と も 高 い コ ス ト の ル ー プ で あ っ た「 literm_k」に つ い て の 性 能 推 定 結 果 をス ラ イ ド p.19 に 示 す 。 FX100 で の 詳 細 PA 実 行 時 間 は 前 述 の と お り 、 128reg、 3

2reg と も に ほ と ん ど 変 化 は な く 、 メ モ リ バ ン ド 幅 ネ ッ ク の ル ー プ で あ っ た 。 一
方 、 性 能 推 定 ツ ー ル に よ る Post-FX100 の 推 定 値 （ ス レ ッ ド 数 は 16 か ら 12 に 変更 ）は 、FX100 で の 結 果（ 約 69 秒 ）か ら 約 28%の 高 速 化 と な る 約 50 秒 で あ っ た 。
メモリビジー時間がほぼ全経過時間を律速しており、メモリバンド幅ネックのループであることにより、演算コストではなくメモリ側のコストで性能が律速されていることが分かる。
一 方 、 も う 一 つ の 高 コ ス ト ル ー プ で あ っ た 「 exb_realspcal」 に 対 す る 前 述 のコ ス ト 分 析 で は 、 レ ジ ス タ 数 の 現 象 に よ る 影 響 は ほ と ん ど な く 、 「 literm_k」と 同 じ く メ モ リ バ ン ド 幅 ネ ッ ク の ル ー プ で あ る こ と が 判 明 し て い た が 、 L1 キ ャ
ッ シ ュ ミ ス 率 が 理 論 値 よ り も 大 き か っ た 。 ス ラ イ ド p.21 に 示 し た Post-FX100で の 性 能 推 定 結 果 で は 、約 20%の 高 速 化 を 示 し て い る が 、や は り 、キ ャ シ ュ の ス
ラ ッ シ ン グ を 反 映 し て 、L2 ビ ジ ー と な っ て い る た め 、デ ー タ を L1 キ ャ ッ シ ュ に適切に乗せるようなチューニングが必要となる。

4.1.5

高コストループのチューニング

GKV カ ー ネ ル コ ー ド に お い て 2 番 目 に 高 い 計 算 コ ス ト を 生 じ て い た 「 exb_rea
lspcal」 に 対 し て 、 Post-FX100 向 け の 最 適 化 に 取 り 組 ん だ 。 前 述 の よ う に 、 この ル ー プ は 、 高 い L1 キ ャ ッ シ ュ ミ ス 率 の た め に L2 ビ ジ ー と な っ て い た た め 、 L
1 キ ャ ッ シ ュ に 載 せ る た め の チ ュ ー ニ ン グ が 必 要 で あ る 。 そ こ で 、ス ラ イ ド p.2
3 に 示 す よ う に 、 該 当 ル ー プ に つ い て 、 （ nx,ny） の 2 次 元 ル ー プ 内 の 演 算 を 2
つ に 分 割 し 、 そ れ ぞ れ を 別 々 の ル ー プ に 分 割 す る こ と で 、 L1 キ ャ ッ シ ュ ミ ス の
低 減 を 期 待 し た 。 結 果 は 、 FX100 の PA 情 報 （ ス ラ イ ド p.24） に 示 す 通 り 、 128r
eg お よ び 32reg コ ン パ イ ラ と も に 、 ル ー プ 分 割 し な い 場 合 に 5%以 上 あ っ た L1D
キ ャ ッ シ ュ ミ ス 率 が 、 2.7%弱 に 低 減 し て お り 、 理 論 値 （ 3.125%） 以 下 に 改 善 され た こ と が 確 か め ら れ た 。こ の PA 情 報 に 基 づ い て 、Post-FX100 で の 性 能 を 推 定
し た 結 果 、 ル ー プ 分 割 し な い 場 合 の FX100 で の 結 果 （ 約 4.7 秒 ） に 対 し て 、 35%の 高 速 化 と な る 約 3.1 秒 ま で 性 能 向 上 す る こ と が 確 か め ら れ た 。詳 細 PA 情 報 から は 、 L2 ア ク セ ス 待 ち が 改 善 さ れ て 、 ス ラ ッ シ ン グ が 解 消 さ れ て い る こ と が 分
か る （ ス ラ イ ド p.26） 。 推 定 に 関 し て も 、 L2 ビ ジ ー の 上 限 に 近 い 性 能 と な っ て
いる。

4.1.6

まとめ

プ ラ ズ マ 乱 流 コ ー ド GKV の カ ー ネ ル コ ー ド を 対 象 と し て 、Post-FX100 に お け る
性 能 の 推 定 を 行 っ た 。GKV は 、プ ラ ズ マ 流 体 解 析 コ ー ド と 異 な り 、一 本 の 磁 力 線に 沿 っ た 局 所 フ ラ ッ ク ス・チ ュ ー ブ 上 で シ ミ ュ レ ー シ ョ ン を 行 う 。5 次 元 位 相 空
間は、空間離散化として磁力線垂直方向の 2 次元空間に対してはフーリエ・ス
ペクトル法、磁力線平行方向の実空間、および 2 次元の速度空間方向について
は 、 4 次 ま た は 5 次 精 度 差 分 法 を 用 い て お り 、 時 間 積 分 に は 4 次 精 度 陽 的 Runge

-Kutta-Gill 法 を 採 用 し て い る 。 本 性 能 推 定 は 、 GKV カ ー ネ ル コ ー ド を 対 象 に １ノ ー ド １ プ ロ セ ス で の 単 体 性 能 を ベ ー ス に 進 め た 。ま ず 、FX100 上 で コ ス ト 分 析を 行 っ た 結 果 、 「 literm_k」 お よ び 「 exb_realspcal」 の サ ブ ル ー チ ン が そ れ ぞれ 全 体 の 約 53.7%、 約 13.6%を 占 め て い た た め 、 こ の 2 つ を 評 価 対 象 に し た 。
対 象 の ル ー プ は 双 方 と も に 、128reg か ら 32reg へ の コ ン パ イ ラ の 変 更 に よ る
性能の劣化はほとんど見られず、メモリビジー率が高く、メモリバンド幅ネッ
ク の ル ー プ で あ る こ と が 分 か っ た 。 「 literm_k」 で の L1D キ ャ ッ シ ュ の ミ ス 率は 理 論 値 付 近 で あ っ た 一 方 、「 exb_realspcal」 で は 、 理 論 値 を 大 き く 超 え て いたため、スラッシングが発生している可能性を確認した。
性 能 予 測 ツ ー ル を 用 い た Post-FX100 で の 性 能 推 定 の 結 果 、 「 literm_k」 は 、
上述のとおり、メモリバンド幅ネックのループであることにより、演算コスト
よ り も メ モ リ 側 の コ ス ト で 性 能 が 律 速 さ れ 、FX100 か ら 約 28%の 高 速 化 が 確 認 され た 。 も う 一 方 の 「 exb_realspcal」 に つ い て は 、 同 様 に メ モ リ バ ン ド 幅 ネ ッ クの ル ー プ で あ る と 同 時 に 、L1 キ ャ ッ シ ュ ミ ス 率 が 理 論 値 よ り も 大 き か っ た た め 、Post-FX100 で の 性 能 推 定 の 結 果 は 、 約 20%の 高 速 化 を 示 し 、 キ ャ シ ュ の ス ラ ッ
シ ン グ に よ り L2 ビ ジ ー と な り 、 そ こ が 性 能 を 律 速 し て い る こ と が 分 か っ た 。
次 に 、 「 exb_realspcal」 に お け る キ ャ シ ュ の ス ラ ッ シ ン グ を 回 避 す る た め 、
ループ内の演算を 2 つに分割し、それぞれを別々のループに分割することによる L1 キ ャ ッ シ ュ ミ ス の 低 減 を 期 待 し た チ ュ ー ニ ン グ を 行 っ た 。チ ュ ー ニ ン グ 後の 性 能 は 、 L1D キ ャ ッ シ ュ ミ ス 率 が 理 論 値 以 下 に 改 善 さ れ 、 Post-FX100 で の 性
能 は 最 適 化 前 に 対 し て 35%の 向 上 を 確 認 し た 。 詳 細 PA 情 報 か ら 、 L2 ア ク セ ス 待ちが改善されて、スラッシングが解消されていることが分かった。
今回の性能推定では、高コストとなっているメモリバンド幅ネックのループ
は ソ フ ト ウ ェ ア ・ パ イ プ ラ イ ン 化 や SIMD 化 が な さ れ て い て 、 性 能 的 な 課 題 が な
い 中 、 高 コ ス ト ル ー プ の 一 つ で あ っ た 「 exb_realspcal」 で は 、 ル ー プ 分 割 を 施
し た 最 適 化 後 の 性 能 で は 、 L2 ビ ジ ー の ル ー プ と な っ た 。 今 後 の さ ら な る 最 適 化としては、ブロッキングにより性能向上の可能性がある。

メニーコア時代のアプリ性能検討WG
LHD turbulencesimulation by GKV

核融合プラズマ乱流コードGKVのPost-FX100性能推定
沼波 政倫 (自然科学研究機構 核融合科学研究所)
片桐 孝洋 (名古屋大学 情報基盤センター)
青木 正樹 (富士通株式会社)1

アウトライン プログラム概要
 ポストFX100向けチューニング用コンパイラによるコスト分析
 Post-FX100における性能の推定
 高コストループのチューニング
 まとめ

2

GKVコード

T.-H. Watanabe, NF (2006), M. Nunami, PFR (2010, 2015),
S. Maeyama, CPC (2013), M. Nakata, CPC (2015)• GyroKinetic Vlasov code : One of major gyrokinetic codes in Japan• Local delta-f model in a flux-tube geometry (shorter time than full-f model).
• Solve gyrokinetic Vlasov-Poisson-Ampere equations- Reduce computations byeliminating gyration.
- CFD calculation in 5D phase space.
- Eulerian (Continuum) solver:spectral in 2D (kx,ky)-space (2D-FFT),
finite-difference in 3D (z, v||, μ)-space.

• High scalability for wide scale
range simulations includingion- and electron-scale.

3

Comm.-Compt. overlap w/ Assistant Core Comm. – Compt. overlap using “OpenMP Master”+MPI
 MPI comm. on OMP Master thread and the independent calculations on
the other threads overlaps.
Idomura, Int’l J.HPC Appl. 2014,
 Utilizing the dynamic scheduling in OMP.
Maeyama, Parallel Compt. 2015 Unable: “Static” scheduling in OMP.

The code optimization and the performance
evaluation are in progress at NIFS.
4AC

16cores MPI comm. on AC and the indep. calculations on“ALL threads” overlaps.
 Enable: fully non-blocking communications,e.g., Send/Recv, AlltoAll, Allreduce.
 Enable: ”Static” schedulingor “chunk-size-optimized dynamic” sched.

16cores Comm. – Compt. overlap using Assistant Core(AC) + MPI on FX100SPARC64TM XIfx

Comm.-Compt. overlap w/ Assistant Core
Master  AC

Integration of calc. costs
w/ Master thread(Dynamic)

thread0w/ AC(Dynamic)

thread0 Almost ideal improvement of ~6.25% = (1/16) by AC.
5

評価環境 FX100 (ハードウェア、ソフトウェアの情報)
FX100 諸元

C
P
U

名称
動作周波数
コア数

SPARC64 XIfx
1.975GHz
32コア＋アシスタントコア

キャッシュメモリ容量

L1I$: 64KB/コア(4way)
L1D$: 64KB/コア(4way)
L2$ : 24MB/CPU

理論ピーク性能
主記憶-CPU間帯域
CPU数、コア数
ー

1011GFlops/CPU(倍精度)2022GFlops/CPU(単精度)(240+240)GB/s/CPU
1CPU/ノード(16コアｘ2CMG＋アシスタントコア)/ノード
1011GFlops/ノード(倍精度)2022GFlops/ノード(単精度)

ノ理論ピーク性能ド

主記憶容量

32GB/ノード

主記憶-CPU間帯域
ノード間通信帯域(240+240)GB/s
100GB/s

OS の名称とバージョン

FX100向けOS 2.1.1

FUJITSU CONFIDENTIAL

6

Copyright 2018 FUJITSU LIMITED

測定条件と分析する高コストループの選定
 FX100 の測定条件
項目

条件

コンパイラ

ポストFX100向けチューニング用コンパイラ

翻訳時オプション

ベース ： -Kfast -Kparallel -Kocl -Cpp -X932reg ： -Knounroll 、ループ分割時 ： -Kloop_nofusion

パラメータ

ループ回転数の指定 (並列化の軸は48回転) ： nm=47 、 nz=24( 前回 16回転時 ： nm=15 、 nz=16 )スレッド数

16スレッド

実行条件

1ノード – 1プロセス (MPI 無し)

 分析する高コストループの選定
 16スレッド(128reg)のプロファイラ情報(Thread 0 の抜粋)から、上位2つのループ(合計で約67%のコスト)を分析する。
Cost% Operation (S) Barrier% Nest Kind Exec Start End---------------------------------------------------------------------------1199 100.0000
16.1980
14 1.1676--- --- ----------------------------------------------------------------------------644
53.71148.70020 0.0000
5 DO
AUTO
62 69
164
13.67812.21560 0.0000
4 DO
AUTO
280 30280
6.67221.08080 0.0000
4 DO
AUTO
71 99
24
2.00170.32420 0.0000
4 DO
AUTO
102 106
FUJITSU CONFIDENTIAL

7

Thread 0
gkv̲advnc.literm̲k.̲PRL̲1̲
①
gkv̲exb.exb̲realspcal.̲PRL̲1̲
②
gkv̲colli.colli̲lb̲model.̲PRL̲1̲
gkv̲exb.exb̲input.̲PRL̲1̲
Copyright 2018 FUJITSU LIMITED

高コストループ①
literm_k の分析
 FX100 の分析
 最適化の状況
 128reg と 32regコンパイラ の比較

FUJITSU CONFIDENTIAL

Copyright 2018 FUJITSU LIMITEDループ①：最適化の状況 (高コスト部のソースリスト)
 高コスト1位：gkv̲advnc.literm̲k.̲PRL̲1̲ (Kernel literm̲k 区間)
50

51

53

59
61

62
63
64
65
66
67
68
69
70
71
72
73
75

!$OMP parallel do private(iv,iz,my,mx)<<< Loop-information Start >>><<< [PARALLELIZATION]<<<
Standard iteration count: 2(ネスト2と3のループはコンパイラでループ交換される)<<< Loop-information End >>>1 pp
do im = 0, nm⇒ 並列化の軸となるループ：48回転 (0〜48) nm は 47<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
INTERCHANGED(nest: 3)(ネスト2と3のループはコンパイラでループ交換される)<<< Loop-information End >>>2 p
do iv = 1, 2*nv
⇒ 16回転 (1〜16) nv は 8<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
INTERCHANGED(nest: 2)<<< Loop-information End >>>3 p
do iz = -nz, nz-1
⇒ 48回転 (-24〜23) nz は 24
4 p
do my = ist̲y, iend̲y⇒ 21回転 (0〜20) ist̲y は 0 , iend̲y は 20<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)<<<
SOFTWARE PIPELINING(ソフトウェアパイプライニングは、ループの回転数が84回以上で実行)<<< Loop-information End >>>(32regコンパイラでは20回以上で実行)
5 p 3vdo mx = -nx, nx⇒ 171回転 (-85〜85) nx は85
5 p 3vlf(mx,my,iz,iv,im) =&5- ui * kvd(mx,my,iz,iv,im) * ff(mx,my,iz,iv,im)&5- cs1 * fmx(iz,iv,im) * (&5+ ui * kvd(mx,my,iz,iv,im) * psi(mx,my,iz,im) &5- ui * kvs(mx,my,iz,iv,im)&5* ( psi(mx,my,iz,im) - cs2 * vl(iv) * chi(mx,my,iz,im)))5 p 3vend do4 p
end do3 p
end do2 p
end do1 p
end do
倍精度複素数 (complex)
倍精度実数 (real)!$OMP end parallel do

FUJITSU CONFIDENTIAL

9

Copyright 2018 FUJITSU LIMITEDループ①：128reg と 32regコンパイラ の比較(1/2)
 PAの実行時間で比較
 128reg と 32regコンパイラで、性能差はほぼ出ていない事から、レジスタ数減少の影響は、小さいと考えられる。
literm_k
FX100 (PA実行時間)

レジスタ数

スレッド数

実行時間(sec)

比率128reg

16

69.13

1.00

32reg

16

69.09

0.99

• 32regでは128regと同様にSWPを動作させるため、-Knounroll を指定

FUJITSU CONFIDENTIAL

10

Copyright 2018 FUJITSU LIMITEDループ①：128reg と 32regコンパイラ の比較(2/2)
 PA情報で比較 (Memory、Cache) 128reg と 32regコンパイラ共にメモリビジー率が高く、メモリバンド幅ネックのループである。
 L1キャッシュのミス率も、理論値(3.15%)付近である。
literm_k

FX100(PA情報)

レジスタ
数

スレッド
数

メモリ
ビジー率

メモリ
スループット(GB/sec)

L1D ミス率(/ロード・
ストア数)

L1D ミス
dm率(/L1D ミス数)

128reg1696%136.68

3.22%

7.31%

32reg

16

96%

136.76

3.18%

7.62%

FX100
128reg
16t

FX100
32reg
16t

浮動小数点ロードL2アクセス待ち

浮動小数点ロード
メモリアクセス待ち

メモリ・キャッシュビジー待ち

FUJITSU CONFIDENTIAL

11

Copyright 2018 FUJITSU LIMITED

高コストループ②
exb_realspcal の分析
 FX100 の分析
 最適化の状況
 128reg と 32regコンパイラ の比較

FUJITSU CONFIDENTIAL

Copyright 2018 FUJITSU LIMITEDループ②：最適化の状況 (高コスト部のソースリスト)
 高コスト2位：gkv̲exb.exb̲realspcal.̲PRL̲1̲ (Kernel exb̲realspcal 区間)

274
277

1
1

p

278

2 pp

279

3

p

280
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
306
307

4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
4
3
2
1
1

p
p

p
p
p
p

<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
INTERCHANGED(nest: 2)(ネスト1と2のループはコンパイラでループ交換される)<<<
PREFETCH
: 2
<<<
vl: 2
<<< Loop-information End >>>do iv = 1, 2*nv
⇒ 16回転 (1〜16) nv は 8!$OMP parallel do private(mx,my)<<< Loop-information Start >>><<< [PARALLELIZATION]<<<
Standard iteration count: 2<<< [OPTIMIZATION]<<<
INTERCHANGED(nest: 1)(ネスト1と2のループはコンパイラでループ交換される)<<< Loop-information End >>>do iz = -nz, nz-1
⇒並列化の軸となるループ： 48回転 (-24〜23) nz は 24<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
PREFETCH: 20<<<
wkexw̲xw: 10, wkdf2̲xw: 10<<< Loop-information End >>>do mx = ist̲xw, iend̲xw⇒ 128回転 (0〜127) ist̲xw は 0 , iend̲xw は 127<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)※ 128reg ではフルアンロールが動作。32regコンパイラでは、<<<FULL UNROLLING※ アンロールではなくSWPL が適用される。(回転数が36回以上で実行)<<< Loop-information End >>>fv
do my = 0, nyw⇒ 65回転 (0〜64) nyw は 64
fv
wkdf1̲xw(my,mx,iz,iv) =&cmplx(
real ( wkdf1̲xw(my,mx,iz,iv), kind=DP ) &* real ( wkeyw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbyw̲xw(my,mx,iz),
kind=DP ) &real ( wkdf2̲xw(my,mx,iz,iv), kind=DP ) &* real ( wkexw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbxw̲xw(my,mx,iz),
kind=DP ) &,aimag ( wkdf1̲xw(my,mx,iz,iv)) &* aimag ( wkeyw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbyw̲xw(my,mx,iz)) &aimag ( wkdf2̲xw(my,mx,iz,iv)) &* aimag ( wkexw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbxw̲xw(my,mx,iz)) &, kind=DP ) * cef
fv
end doend doend do!$OMP end parallel doend do

FUJITSU CONFIDENTIAL13

倍精度複素数 (complex)
倍精度実数 (real)

Copyright 2018 FUJITSU LIMITEDループ②：128reg と 32regコンパイラ の比較(1/2)
 PAの実行時間で比較
 128reg と 32regコンパイラで、性能差は大きく出ていない事から、レジスタ数減少の影響は、小さいと考えられる。
exb̲realspcal
FX100 (PA実行時間)

レジスタ数

スレッド数

実行時間(sec)

比率128reg

16

4.72

1.00

32reg

16

4.57

0.96

• 32regでは、SWPを動作させるため、-Knounroll を指定
128regは、フルアンローリング動作

FUJITSU CONFIDENTIAL

14

レジスタ数で最適化が異なる
128reg：SIMD + FULL UNROLLING
32reg：SIMD + SWPL

Copyright 2018 FUJITSU LIMITEDループ②：128reg と 32regコンパイラ の比較(2/2)
 PA情報で比較 (Memory、Cache) 128reg と 32regコンパイラ共にメモリビジー率が高く、メモリバンド
幅ネックのループであるが、L1D ミス率が理論値(3.15%)よりも高い
ため、スラッシングが発生していると考えられる。
exb̲realspcal

FX100(PA情報)

レジスタ数

スレッド
数

メモリ
ビジー率

メモリ
スループッ
ト(GB/sec)

L1D ミス率(/ロード・
ストア数)

L1D ミス
dm率(/L1D ミス数)

128reg16

89%

119.745.12%

48.80%

32reg

16

85%

115.82

5.24%

50.67%

 スラッシングの対応として、ループ分割によるチューニングを実施する。 (前回発表時と同様)

FX100
128reg
16t

FX100
32reg
16t

浮動小数点ロードL2アクセス待ち

メモリ・キャッシュ
ビジー待ち

FUJITSU CONFIDENTIAL

15

Copyright 2018 FUJITSU LIMITEDGKV
Post-FX100推定
 ポストFX100のアーキテクチャ
 性能予測ツールの概要
 高コストループ① literm_k の分析
 高コストループ② exb_realspcal の分析

FUJITSU CONFIDENTIAL

16

Copyright 2018 FUJITSU LIMITED

性能予測ツールの概要
アプリコード（カーネル部）
FX100：32regコンパイラを使用し詳細プロファイラ実行でＰＡ情報を取得

PA情報：詳細プロファイラの出力結果(csvファイル）
性能予測ツール
性能予測値
 性能予測ツール概略
 FX100のＰＡ情報を入力とする
 ハード、命令セット／コンパイラによる性能向上効果を算出する
 性能予測値を出力する

FUJITSU CONFIDENTIAL

17

Copyright 2018 FUJITSU LIMITED

高コストループ①
literm_k の分析
 Post-FX100推定
 FX100 の実行時間と Post-FX100推定値の比較
 推定性能について

FUJITSU CONFIDENTIAL

18

Copyright 2018 FUJITSU LIMITEDループ①：Post-FX100推定
 FX100 の実行時間と Post-FX100推定値の比較
literm_k
FX100 (PA実行時間)
Post-FX100 (推定値)

レジスタ数

スレッド数

実行時間(sec)

比率128reg

16

69.13

1.00

32reg

16

69.09

0.99

32reg

12

50.46

0.72

FX100
128reg
16t

FX100
32reg
16t

P-FX100
32reg 12t推定

 推定性能について メモリバンド幅ネックのループであり、演算のコストよりも、メモリ側のコストで性能が律速するため、Post-FX100への移行に関して、特に問題はない。
ビジー率
Post-FX100

メモリビジー メモリRビ メモリWビ
L2ビジー率 L1ビジー率
率
ジー率
ジー率

91%

70%

37%

100%

浮動小数点ロードL2アクセス待ち57%

Post-FX100推定 ： 32reg 12t asis
浮動小数点ロード
メモリアクセス待ち

メモリ・キャッシュ
ビジー待ち

ビジー待ち

FUJITSU CONFIDENTIAL

19

Copyright 2018 FUJITSU LIMITED

高コストループ②
exb_realspcal の分析
 Post-FX100推定
 FX100 の実行時間と Post-FX100推定値の比較
 推定性能について

FUJITSU CONFIDENTIAL

20

Copyright 2018 FUJITSU LIMITEDループ②：Post-FX100推定
 FX100 の実行時間と Post-FX100推定値の比較
exb̲realspcal
FX100 (PA実行時間)
Post-FX100 (推定値)

レジスタ数

スレッド数

実行時間(sec)

比率128reg

16

4.72

1.00

32reg

16

4.57

0.96

32reg

12

3.78

0.80

 推定性能について

レジスタ数で最適化が異なる
128reg：SIMD+FULL UNROLLING
32reg：SIMD+SWPL
FX100
128reg
16t

FX100
32reg
16t

P-FX100
32reg 12t

推定

 メモリバンド幅ネックのループであるが、キャシュのスラッシングにより、 Post-FX100の推定ではL2ビジーとなるため、L1キャッシュに乗るようなチューニング(ループ分割など)を実施すると、性能改善が期待できる。 (前回発表時と同様)
ビジー率

メモリビ
ジー率

Post-FX100

47%

メモリRビ メモリWビ L2ビジー L1ビジー
ジー率
ジー率
率
率

34%

24%

94%

72%

Post-FX100推定 ： 32reg 12t asis

浮動小数点ロードL2アクセス待ち

メモリ・キャッシュ
ビジー待ち

FUJITSU CONFIDENTIAL

21

Copyright 2018 FUJITSU LIMITED

高コストループ②
exb_realspcal のチューニング
 チューニング(ループ分割)の実施
 最適化の状況
 Post-FX100推定

FUJITSU CONFIDENTIAL

Copyright 2018 FUJITSU LIMITEDループ②：ループ分割：最適化の状況
 高コスト2位を手動でループ分割
266
268
269
271
272

1 pp
1
2 p
2
3 p

273
274
275
276
277
278
279
280
281
282
283
284
285

4
4
4
4
4
4
4
4
4
4
3
2
3

p
p

v
v

p
p

v

286
287
288
289
290
291
292
293
294
295
296
297
299
300

4
4
4
4
4
4
4
4
4
4
3
2
2
1

p
p

v
v

p
p

v

p

p

<<< Loop-information Start >>><<< [PARALLELIZATION]<<<
Standard iteration count: 2<<< Loop-information End >>>do iz = -nz, nz-1
! 手動でループ交換 ※コンパイラで自動的に交換された場合、自動並列化の回転軸が変化するため。
!$OMP parallel do private(mx,my)do iv = 1, 2*nv
! 手動でループ交換!-------------------------------------------------------------------do mx = ist̲xw, iend̲xw ! 手動でループ分割<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)<<<
SOFTWARE PIPELINING! 32regではレジスタが不足するため -Knounroll を指定し、最適化を FULL UNROLLING から SWPL に変更<<< Loop-information End >>>do my = 0, nyw
wkdf1̲xw(my,mx,iz,iv) =&cmplx(
real ( wkdf1̲xw(my,mx,iz,iv), kind=DP ) &* real ( wkeyw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbyw̲xw(my,mx,iz),
kind=DP ) &,aimag ( wkdf1̲xw(my,mx,iz,iv)) &* aimag ( wkeyw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbyw̲xw(my,mx,iz)) &, kind=DP )end doend do!-------------------------------------------------------------------do mx = ist̲xw, iend̲xw ! 手動でループ分割<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)<<<
SOFTWARE PIPELINING! 32regではレジスタが不足するため -Knounroll を指定し、最適化を FULL UNROLLING から SWPL に変更<<< Loop-information End >>>do my = 0, nyw
wkdf1̲xw(my,mx,iz,iv) = (wkdf1̲xw(my,mx,iz,iv) +&cmplx(real ( wkdf2̲xw(my,mx,iz,iv), kind=DP ) &* real ( wkexw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbxw̲xw(my,mx,iz),
kind=DP ) &,aimag ( wkdf2̲xw(my,mx,iz,iv)) &* aimag ( wkexw̲xw(my,mx,iz) - cs1 * vl(iv) &* wkbxw̲xw(my,mx,iz)) &, kind=DP )) * cef
end do
倍精度複素数 (complex)
end do
倍精度実数 (real)!-------------------------------------------------------------------end do!$OMP end parallel do

FUJITSU CONFIDENTIAL

23

Copyright 2018 FUJITSU LIMITEDループ②：ループ分割：Post-FX100推定(1/2)
 ループ分割に関する評価
 FX100のPA情報の比較 (Memory、Cache) ループ分割によりキャシュのスラッシングが解消し、L1Dミス率が理論値(3.125%)以下に改善された。
exb_realspcal

レジスタ数

128reg

スレッド数

16

FX100
(PA情報)
32reg

16

ループ分割

L1D ミス率(/ロード・ストア数)

L1D ミスdm率(/L1D ミス数)

なし

5.12%

48.80%

あり

2.69%

17.76%なし

5.24%

50.67%

あり

2.65%

14.04%

 L1Dキャッシュスラッシングの目安
L1D ミス率(/ロード・ストア数)
単精度 1.5625%以上
倍精度 3.125%以上

FUJITSU CONFIDENTIAL

単精度:1/64（64回に1回ミスをする）
倍精度:1/32（32回に1回ミスをする）

24

Copyright 2018 FUJITSU LIMITEDループ②：ループ分割：Post-FX100推定(1/2)
 ループ分割に関する評価
 Post-FX100推定の性能予測(推定時間)の比較
 ループ分割により、推定性能が改善された。
exb_realspcal

レジスタ数

スレッド数

128reg16

FX100(PA実行時間)
32reg
Post-FX100(推定値)

FUJITSU CONFIDENTIAL

32reg

16

12

ループ分割

推定時間(sec)

比率

なし

4.72

1.00

あり

4.91

1.03

なし

4.57

0.96

あり

4.85

1.02

なし

3.78

0.80

あり

3.10

0.65

25

Copyright 2018 FUJITSU LIMITEDループ②：ループ分割：Post-FX100推定(2/2)
 ループ分割時のPost-FX100推定PA情報の比較
 ループ分割により、スラッシングが解消され、L2アクセス待ち が改善した。

ループ2分割 推定に関しても、L2ビジーの上限に近い性能となっている。
FX100
128reg16t

Post-FX100推定 ： 32reg 12t ループ2分割

ビジー率

FX100
32reg
16t

P-FX100
32reg 12t

推定

メモリビ メモリRビ メモリWビ L2ビジー L1ビジー
ジー率 ジー率 ジー率率率

Post-FX100 65%

49%

29%

94%

72%

浮動小数点ロード
L1Dアクセス待ち

浮動小数点ロードL2アクセス待ち

浮動小数点ロード
メモリアクセス待ち

メモリ・キャッシュビジー待ち

FUJITSU CONFIDENTIAL

26

Copyright 2018 FUJITSU LIMITED

まとめ プラズマ乱流コードGKVのPost-FX100における性能の推定を行った ポストFX100向けコンパイラによるコスト分析
 「literm_k」、「exb_realspcal」が全体の約67%を占める
 レジスタ数の減少(128reg→32reg)による性能劣化はない
 メモリビジー率が高く、基本的にメモリバンド幅ネック
 「exb_realspcal」では、スラッシングが発生

 性能推定ツールによるPost-FX100における性能推定
 「literm_k」は、FX100から約28%の高速化となり、メモリ側コストで性能が律速
 「exb_realspcal」はキャシュ・スラッシングにより約20%の高速化にとどまる

 高コストループのチューニング
 ループを分割し、L1キャッシュミスの低減を期待したチューニングを実施
 最適化後、約35%の性能向上を確認
27
 L2アクセス待ちが改善されて、スラッシングが解消

4.2 宇宙プラズマ５次元ブラソフコード Vlasov5 の性能測定および推定

4.2.1

名古屋大学宇宙地球環境研究所

梅田 隆行

富士通株式会社

内藤 俊也

はじめに

希薄で無衝突状態にある宇宙プラズマは、磁気流体力学、イオン運動論、電子 運 動 論 な ど の 様 々 な 時 空 間 ス ケ ー ル で 物 理 現 象 を 記 述 で き る 。 Vlasov コ ー ド
は 、宇 宙 プ ラ ズ マ の 第 一 原 理 で あ る Vlasov（ 無 衝 突 Boltzmann）方 程 式 と Maxwell
方程式により、プラズマ中の全ての電磁・静電波動と荷電粒子の運動との相互
作 用 を 解 く シ ミ ュ レ ー シ ョ ン 手 法 で あ る (文 末 資 料 p.1、 以 降 ペ ー ジ 番 号 の み 表記 す る )。Vlasov 方 程 式 は 、位 置 及 び 速 度 の 関 数 と し て 表 さ れ る 位 相 空 間 分 布 関
数の保存則である。現実の世界では、位置及び速度は共に 3 次元であるが、既存のコンピュータシステムにおいて 6 次元問題を扱うのは困難であるため、本プ ロ ジ ェ ク ト で 用 い る Vlasov5 で は 実 空 間 （ 位 置 ） を 2 次 元 、 速 度 空 間 を 3 次
元 と し た 5 次 元 問 題 を 扱 う 。 本 稿 で は 、 前 マ ル チ コ ア ク ラ ス タ 性 能 WG 以 降 で 行っ た プ ロ グ ラ ム の 改 良 点 に つ い て ま と め 、 さ ら に ポ ス ト FX100 で の 性 能 を 推 定した。
4.2.2 プ ロ グ ラ ム 概 要
Vlasov コ ー ド は 、 Maxwell 方 程 式 に よ る 電 磁 場 の 更 新 と Vlasov 方 程 式 に よ る
分 布 関 数 の 更 新 か ら 成 る 。Maxwell 方 程 式 は 電 磁 場 の 解 析 で 広 く 用 い ら れ て い る
FDTD（ Finite Difference Time Domain）法 に よ り 、そ の 時 間 発 展 を 解 き 進 め る 。
Maxwell 方 程 式 の 計 算 負 荷 は Vlasov 方 程 式 の 計 算 負 荷 に 対 し て 通 常 0.1%未 満 であるため、今回は性能評価の対象外とした。
Vlasov 方 程 式 に よ る 分 布 関 数 の 更 新 は 、演 算 子 分 離 法 に よ り 以 下 の よ う に 実 空
間（位置方向）の移流、電場による速度空間の移流及び磁場による速度空間の回 転 の ３ つ の 部 分 に 分 解 さ れ 、 そ れ ぞ れ の 式 が ① 実 空 間 の 移 流 （ position カ ーネ ル ） 、 速 度 空 間 の ② 移 流 （ velocity_e カ ー ネ ル ） 及 び ③ 回 転 （ velocity_b カー ネ ル ） に 対 応 す る (p.2)。
∂𝑓𝑠
∂𝑡

∂𝑓𝑠
∂𝑡

∂𝑓

+ 𝐯 ∙ ∂𝐫𝑠 = 0
+

∂𝑓𝑠
∂𝑡

𝑞𝑠
𝑚𝑠

+

𝐄∙

𝑞𝑠
𝑚𝑠

∂𝑓𝑠
∂𝐯

(1)
(2)

=0

(𝐯 × 𝐁) ∙

∂𝑓𝑠
∂𝐯

=0

(3)

式 (1)及 び (2)は 線 形 移 流 方 程 式 で あ り 、演 算 子“ 非 ”分 離 型 の 保 存 型 解 法 を 用
い る 。 ま た 式 (3)は 回 転 流 方 程 式 で あ り 、 back-substitution 法 と 呼 ば れ る 演 算
子分離型の解法を用いる。またこれらの方程式を解く上で、物理量の保存、解の無振動性、正値性の保証を満たす独自の 5 次精度保存型無振動スキームを用い て い る (p.3)。

速 度 空 間 の 移 流 （ velocity_e カ ー ネ ル ） の 概 要 を プ ロ グ ラ ム 1 に 示 す (p.4)。
実 空 間 x,y 方 向 の 添 え 字 を i,j、 速 度 空 間 vx,vy,vz の 添 え 字 を l,m,n と す る 。
こ の プ ロ グ ラ ム で は 、 vx,vy,vz の 各 方 向 の 1 次 元 数 値 流 束 を 計 算 し た の ち 3 次
元 数 値 流 束 を 合 成 し て い る (p.4)。分 布 関 数 デ ー タ の 定 義 は 、モ ー メ ン ト 計 算 (速度 空 間 の 積 分 )を 高 速 に 行 う た め に f(vx,vy,vz,x,y)と な っ て お り 、 velocity_e
カ ー ネ ル で は 数 値 流 束 デ ー タ が l,n,m に 依 存 す る た め に 3 次 元 配 列 と な る 。 一
方 、 実 空 間 の 移 流 （ position カ ー ネ ル ） で は 数 値 流 束 デ ー タ が i,j に 依 存 す る
た め に (l,m,n,i,j)の 5 次 元 配 列 と な る べ き で あ る が 、分 布 関 数 デ ー タ を 転 置 して f(x,y,vx,vy,vz)と 定 義 す る こ と に よ り 、 2 次 元 配 列 と な る 。 ま た こ の 配 列 の転 置 に よ り 、position カ ー ネ ル は プ ロ グ ラ ム 2 の velocity_e カ ー ネ ル と ほ と ん
ど 同 じ と な る (p.6)。
プ ロ グ ラ ム 1： velocity_e カ ー ネ ル の 概 要12
3

DO j=1,Ny
DO i=1,Nx
DO n=0,Nvz

4

DO m=0,Nvy

5

DO l=0,Nvx

6dfx(l,m,n)=... !vx 方 向 の フ ラ ッ ク ス の 計 算

7dfy(l,m,n)=... !vy 方 向 の フ ラ ッ ク ス の 計 算

8dfz(l,m,n)=... !vz 方 向 の フ ラ ッ ク ス の 計 算

9END DO

10

END DO

11

END DO

12

DO n=1,Nvz

13

DO m=1,Nvy

14

DO l=1,Nvx

15f(l,m,n,i,j)=f(l,m,n,i,j) -dfx(l,m,n)+dfx(l-1,m,n) &16-dfy(l,m,n)+dfy(l,m-1,n) &17-dfz(l,m,n)+dfz(l,m,n-1)

18

END DO

19
20
21

END DOEND DOEND DO

22 ND DO

4.2.3 測 定 環 境 お よ び プ ロ グ ラ ム
性 能 測 定 に は 、 名 古 屋 大 学 の FX100 を 使 用 し た 。 な お 、 本 測 定 は 単 一 ノ ー ド
で 行 い 、 2 プ ロ セ ス ・ 16 ス レ ッ ド を 利 用 し た 。 ま た ノ ー ド 間 の 並 列 性 能 の 測 定は行わない。コンパイラオプションは以下のとおりである。
 -Kfast,openmp,ocl,simd=2 –x250
測 定 に は ２ つ の プ ロ グ ラ ム を 使 用 し た 。 プ ロ グ ラ ム ２ は WG 開 始 時 点 で の asis コ ー ド で あ り 、 プ ロ グ ラ ム 1 の 3-11 行 目 の vx,vy,vz の 各 速 度 方 向 の 数 値 流
束 を 計 算 す る プ ロ グ ラ ム で あ る 。 変 数 lv,l0,mv,m0,nv,n0,sx,sy,sz は 、 実 空 間の 変 数 で あ る 電 場 ベ ク ト ル の 符 号 に 依 存 し て +1 ま た は -1 の ど ち ら か の 値 を と り 、
1 行 目 の ル ー プ の 外 で 計 算 さ れ る 。 ま た inv3 は parameter 定 数 1/3 で あ る 。 プ
ロ グ ラ ム ３ は 、 プ ロ グ ラ ム ２ の 30—47 行 目 の 累 積 和 の 計 算 を 複 数 の 1 次 元 ル ープに分割したものである。
4.2.4

累積計算のループ分割

プ ロ グ ラ ム ２ で は 、35—46 行 目 に あ る 3 次 元 配 列 へ の 累 計 計 算 に お い て 配 列 要
素 へ の ア ク セ ス に 依 存 関 係 が あ り 、 30 行 目 の ル ー プ に 対 す る コ ン パ イ ラ の 最 適
化 を 阻 害 し て い た (p.7)。 プ ロ グ ラ ム ３ で は 、 36,41,46,51 行 目 の ル ー プ 内 に お
い て dfx,dfy,dfz の 各 配 列 に ア ク セ ス す る の は そ れ ぞ れ 1 回 で あ り 、 配 列 要 素へのアクセスに対して依存関係が無いため、各ループについてコンパイラの最適 化 が 促 進 さ れ る (p.8)。 た だ し 、 30 行 目 の ル ー プ に お け る gfxy,gfyz,gfzx,gf
xyz の 計 算 結 果 を 、分 割 し た 36,41,46,51 行 目 の ル ー プ に 渡 す 必 要 が あ る 。こ のため、プログラム３ではこれらの変数を 1 次元配列として定義する必要が生じた。
プ ロ グ ラ ム ２ を 2 プ ロ セ ス ・ プ ロ セ ス あ た り の 16 ス レ ッ ド で 実 行 す る と 、 1
タ イ ム ス テ ッ プ あ た り の 計 算 に 約 3.19 秒 掛 か っ た 。 こ の 最 適 化 に よ っ て 、 プ ロ
グ ラ ム ３ で は 計 算 時 間 が 約 2.48 秒 に 短 縮 さ れ 、約 1.23 倍 の 高 速 化 に 成 功 し た 。
ループ分割には配列が増えることによって計算性能が劣化するデメリットもあるが、本例の場合にはコンパイラの最適化による性能改善のメリットの方が大
き か っ た (p.9,10)。
4.2.5 128reg/32reg コ ン パ イ ラ 比 較 お よ び 、 ツ ー ル に よ る ポ ス ト FX100 の 性
能推定
FX100 に 搭 載 さ れ て い る SPARC64XIfx プ ロ セ ッ サ の レ ジ ス タ 数 は 128 で あ る が 、
ポ ス ト FX100 に 採 用 さ れ る A64FX プ ロ セ ッ サ で は レ ジ ス タ 数 が 32 に 削 減 さ れ る 。
本 節 で は 、FX100 上 に 用 意 さ れ た 、レ ジ ス タ 数 の 削 減 を 考 慮 し た 32 レ ジ ス タ 版
富 士 通 コ ン パ イ ラ（ 以 下 、32reg コ ン パ イ ラ と よ ぶ ）を 用 い て プ ロ グ ラ ム ３ を 測
定 し た 。 さ ら に 、 推 定 ツ ー ル に よ り 、 プ ロ グ ラ ム ３ の ポ ス ト FX100 上 で の 性 能を 推 定 し た (p.11)。

32reg コ ン パ イ ラ を 用 い て 翻 訳 し た プ ロ グ ラ ム は 、128reg コ ン パ イ ラ

を 用 い て 翻 訳 し た プ ロ グ ラ ム に 対 し て 、プ ロ グ ラ ム 全 体 で 速 度 が 約 68%に 低下 し た (p.12)。 な お 、 128reg コ ン パ イ ラ を 用 い た 場 合 に は す べ て の カ ー ネ
ル に お い て ocl に よ る 手 動 ル ー プ ア ン ロ ー リ ン グ を 用 い た 時 が 最 も 高 速 であ っ た が 、 32reg コ ン パ イ ラ を 用 い た 場 合 は velocity_e カ ー ネ ル に つ い ての み 大 幅 な 速 度 低 下 が 生 じ た た め 、ocl に よ る 手 動 ル ー プ ア ン ロ ー リ ン グ を外 し た 。な お 、こ の ocl に よ る 手 動 ル ー プ ア ン ロ ー リ ン グ に つ い て は 、ポ ス
ト ペ タ ア プ リ 性 能 WG 成 果 報 告 書 を 参 照 さ れ た い 。 負 荷 の 高 い ル ー プ お よ び
行 は 、関 数 pic5 を 呼 び 出 し て い た (p.13,14,15)。ま た 関 数 pic5 内 で は 組 み
込 み 関 数 を 呼 び 出 す 行 の 負 荷 が 高 い 可 能 性 が あ る こ と が 分 か っ た (p.16)。性
能 劣 化 が 一 番 激 し か っ た velocity_b カ ー ネ ル で は 、 浮 動 小 数 点 演 算 待 ち が増 加 し て お り 、レ ジ ス タ 数 が 少 な く な る こ と で 演 算 待 ち が 生 じ た 可 能 性 が ある (p.17)。
32reg コ ン パ イ ラ に よ る 測 定 結 果 を 性 能 推 定 ツ ー ル に 入 力 す る こ と に よ り 、
ポ ス ト FX100 に お け る コ ー ド の 性 能 を 推 定 す る こ と が で き る 。こ の ツ ー ル に
よ る と 、32reg コ ン パ イ ラ に よ る 測 定 結 果 に 対 し て 、ポ ス ト FX100 で の 実 行
時 間 は 約 1.28 倍 に な る と 推 定 さ れ た 。 た だ し 、 性 能 推 定 ツ ー ル で は ノ ー ド
あ た り 2 プ ロ セ ス ・ 12 ス レ ッ ド で あ る の に 対 し 、 ポ ス ト FX100 に 採 用 さ れる A64FX プ ロ セ ッ サ は 48 コ ア あ る た め 、実 行 時 間 は 約 0.64 倍 に な る と 推 定
さ れ る (p.18)。本 プ ロ グ ラ ム は 演 算 ネ ッ ク で あ り 、性 能 推 定 ツ ー ル が 実 際 の性 能 よ り 低 め の 推 定 値 を 出 す 傾 向 に あ る 。 ル ー プ 内 で ス カ ラ 関 数 pic5 を 呼
ぶ プ ロ グ ラ ム 構 造 か ら 、ベ ク ト ル (配 列 )関 数 を 呼 ぶ プ ロ グ ラ ム 構 造 に 変 更 し 、
関 数 内 で ル ー プ 分 割 を 行 う こ と に よ り (p.19)、約 1.7 倍 の 高 速 化 が 見 込 め る 。
こ れ は 、128reg コ ン パ イ ラ に よ る 測 定 結 果 に 対 し て 、ポ ス ト FX100 の 48 コ
ア で の 実 行 時 間 は 約 0.55 倍 に な る こ と を 意 味 す る 。

4.2.6

まとめ

性 能 チ ュ ー ニ ン グ に よ っ て 、 式 (2)を 解 く プ ロ グ ラ ム 全 体 に お い て 、 約 2.37
倍 の 高 速 化 が な さ れ た 。 以 下 は 、 本 WG の 活 動 で 得 ら れ た 知 見 で あ る 。
 配列に対する累積計算は、ループ分割により配列に対するアクセスの依存性を無くすことができ、ループの最適化及び計算の高速化が可能である。
 レジスタ数が 128 から 32 に削減されることにより、プ ロ グ ラ ム 全 体 で 速 度 が 約68%に 低下する可能性がある。これにより、A64FX プ ロ セ ッ サ の 48 コ ア で は 、
SPARC64 XIfx の 32 コ ア に 対 し て プ ロ グ ラ ム 全 体 の 実 行 時 間 が 約 0.94 倍(=0.64/0.68)に な る と 推 定 さ れ る 。
 翻訳指示文(OCL)によるループアンローリング数の調整や、ソフトウェアパイプラ
イニングの抑制により、性能が改善する可能性がある。

 ループ内でスカラ関数を呼ばずに、サブルーチンによりループそのものをベクトル演算することにより、ループ分割などの最適化手法を適用することが、今後の課題として挙げられる。これにより、約 1.7 倍 の 高 速 化 が 見 込 め る 。

プログラム２
1 do n=0,nvz+1
2
do m=0,nvy+1
3
do l=0,nvx+1
4ffi(l)=1.0d0/ff(l,m,n,i,j)
5
end do6
do l=0,nvx+1
7
hm2=ff(l-lv-lv,m,n,i,j)8
hm1=ff(l-lv ,m,n,i,j)9
hp0=ff(l,m,n,i,j)10
hp1=ff(l+lv ,m,n,i,j)11
hp2=ff(l+lv+lv,m,n,i,j)12
gfx(l)=pic5(hp2,hp1,hp0,hm1,hm2,ex(i,j))13
end do14
do l=0,nvx+1
15
hm2=ff(l,m-mv-mv,n,i,j)
16
hm1=ff(l,m-mv ,n,i,j)
17
hp0=ff(l,m,n,i,j)18
hp1=ff(l,m+mv ,n,i,j)
19
hp2=ff(l,m+mv+mv,n,i,j)20
gfy(l)=pic5(hp2,hp1,hp0,hm1,hm2,ey(i,j))21
end do22
do l=0,nvx+1
23
hm2=ff(l,m,n-nv-nv,i,j)24
hm1=ff(l,m,n-nv ,i,j)25
hp0=ff(l,m,n,i,j)26
hp1=ff(l,m,n+nv ,i,j)
27
hp2=ff(l,m,n+nv+nv,i,j)28
gfz(l)=pic5(hp2,hp1,hp0,hm1,hm2,ez(i,j))29
end do30
do l=0,nvx+1
31
gfxy = abs((gfx(l)*ffi(l))* gfy(l))*0.5d032
gfyz = abs((gfy(l)*ffi(l))* gfz(l))*0.5d033
gfzx = abs((gfz(l)*ffi(l))* gfx(l))*0.5d034
gfxyz = abs((gfx(l)*ffi(l))*(gfy(l)*ffi(l))*gfz(l))*inv3
35dfx(l+l0,m ,n ) = dfx(l+l0,m,n ) +(gfx(l)+(gfxyz-gfxy-gfzx)*sx)
36dfx(l+l0,m ,n+nv) = dfx(l+l0,m,n+nv)-(gfxyz-gfzx)*sx
37
dfx(l+l0,m+mv,n ) = dfx(l+l0,m+mv,n )-(gfxyz-gfxy)*sx
38
dfx(l+l0,m+mv,n+nv) = dfx(l+l0,m+mv,n+nv)+ gfxyz*sx
39
dfy(l ,m+m0,n ) = dfy(l ,m+m0,n ) +(gfy(l)+(gfxyz-gfyz-gfxy)*sy)
40
dfy(l+lv,m+m0,n ) = dfy(l+lv,m+m0,n )-(gfxyz-gfxy)*sy
41dfy(l ,m+m0,n+nv) = dfy(l ,m+m0,n+nv)-(gfxyz-gfyz)*sy42
dfy(l+lv,m+m0,n+nv) = dfy(l+lv,m+m0,n+nv)+ gfxyz*sy
43dfz(l ,m ,n+n0) = dfz(l ,m ,n+n0) +(gfz(l)+(gfxyz-gfzx-gfyz)*sz)
44dfz(l ,m+mv,n+n0) = dfz(l ,m+mv,n+n0)-(gfxyz-gfyz)*sz
45dfz(l+lv,m ,n+n0) = dfz(l+lv,m,n+n0)-(gfxyz-gfzx)*sz
46
dfz(l+lv,m+mv,n+n0) = dfz(l+lv,m+mv,n+n0)+ gfxyz*sz
47
end do48 end do
49 end do

プログラム３
30 do l=0,nvx+1
31
gfxy(l) = abs((gfx(l)*ffi(l))* gfy(l))*0.5d032
gfyz(l) = abs((gfy(l)*ffi(l))* gfz(l))*0.5d033
gfzx(l) = abs((gfz(l)*ffi(l))* gfx(l))*0.5d034
gfxyz(l) = abs((gfx(l)*ffi(l))*(gfy(l)*ffi(l))*gfz(l))*inv3
35 end do
36 do l=0,nvx+1
37dfx(l+l0,m,n)dfx(l+l0,m ,n )+gfx(l)+(gfxyz(l)-gfxy(l)-gfzx(l))*sx
38dfy(l,m+m0,n)dfy(l ,m+m0,n )+gfy(l)+(gfxyz(l)-gfyz(l)-gfxy(l))*sy
39
dfz(l
,m
,n+n0)
dfz(l ,m ,n+n0)+gfz(l)+(gfxyz(l)-gfzx(l)-gfyz(l))*sz
40 end do
41 do l=0,nvx+1
42dfx(l+l0,m ,n+nv) = dfx(l+l0,m ,n+nv)-(gfxyz(l)-gfzx(l))*sx
43dfy(l+lv,m+m0,n ) = dfy(l+lv,m+m0,n )-(gfxyz(l)-gfxy(l))*sy44dfz(l ,m+mv,n+n0) = dfz(l ,m+mv,n+n0)-(gfxyz(l)-gfyz(l))*sz
45 end do
46 do l=0,nvx+1
47
dfx(l+l0,m+mv,n ) = dfx(l+l0,m+mv,n )-(gfxyz(l)-gfxy(l))*sx
48dfy(l ,m+m0,n+nv) = dfy(l ,m+m0,n+nv)-(gfxyz(l)-gfyz(l))*sy49dfz(l+lv,m ,n+n0) = dfz(l+lv,m ,n+n0)-(gfxyz(l)-gfzx(l))*sz
50 end do
51 do l=0,nvx+1
52
dfx(l+l0,m+mv,n+nv) = dfx(l+l0,m+mv,n+nv)+ gfxyz(l)*sx
53dfy(l+lv,m+m0,n+nv) = dfy(l+lv,m+m0,n+nv)+ gfxyz(l)*sy
54
dfz(l+lv,m+mv,n+n0) = dfz(l+lv,m+mv,n+n0)+ gfxyz(l)*sz
55 end do===

運動論方程式：流体方程式の元となる第一原理
粒子間の衝突が無視できるプラズマ（電離気体）と電磁場を扱う
 粒子(Particle-In-Cell)コード

 分布関数(Vlasov)コード
f
f q
f
v
 E  v  B   0
t
x m
v

r = (x,y,z) と v = (vx,vy,vz) のN体
x(n), y(n), z(n), vx(n), vy(n), vz(n)と電磁場 E(i,j,k), B(i,j,k) を扱う

分布関数 f(i,j,k,l,m,n)と電磁場 E(i,j,k), B(i,j,k) を扱う

yx

オイラー（流体）型

電磁場：Maxwell方程式

電磁場格子の中を、粒子が動き回る
ノード間並列化（ロードバランス）に課題

電磁場格子上に分布関数が定義されている高いスケーラビリティ
1コードの概要：Vlasov5D
解いている方程式：
 Maxwell方程式（電磁場の変化と伝搬）
B
1 E
負荷は0.1%未満
ÑE Ñ  B  0J  2c tt Vlasov（無衝突Boltzmann）方程式（プラズマの運動）
f(x, y, z, vx, vy, vz)
f s
f s q s
f s
t

v

x



ms

E  v  B  

v

0

3つに分解
f s
f s
v
0
t
x

f s q s
fE s 0t m sv（位置更新）
position（速度更新：移流）velocity_e

f s
qs
f s
v  B  

0
t
ms
v

（速度更新：回転）velocity_b2

コードの構造 まとめ
移流計算：velocity_e, position
 velocity_e：速度微分、演算が速度空間依存
 position：位置微分、演算が実空間依存
配列の転置(速度,位置)⇒(位置,速度)により、
velocity_eと同じ数値解法(Umeda+ CPC 2009)を適用

回転計算：velocity_b
 速度微分、演算が速度空間依存 剛体回転のため、特殊な演算子分離法(back-substitution法:Schmitz&Grauer CPC 2006)を使用
 vz(loop4)⇒vy(loop5)⇒vx(loop6)⇒vz(loop7)⇒vy(loop8)⇒vz(loop9)の順に計算。4と9、5と8、6と7は同じプログラム構造を持つ

共通解法：5次精度の、保存型・無振動・正値スキーム(Umeda EPS 2008; Umed+ CPC 2012)3f s qs
f
E s  0t ms
v

velocity_e：電場による加速を計算

vz方向流束計算

5次・保存型・無振動・正値スキーム
により1次元流束を計算
重い演算（割り算・流束計算）を1次元ループ分割
ループ間のデータ共有を1次元配列で行う

1次元流束から3次元流束を合成

4

Uy0

Ux0

1次元流束計算のみ：
斜め方向への物理量の輸送は不可⇒多段時間積分(R-K法etc…)
Ux2
Uy1

Uy2

Ux1

1次元流束から多次元流束を合成し、周辺格子点の情報を累積
斜め方向への物理量の輸送を可能に[Umeda et al. CPC 2009]
Ux0=Ux1+Ux2
Uy0=Uy1+Uy2

5

“position” subroutine
𝜕𝑓௦
𝜕𝑓௦
+𝒗ȉ
=0
𝜕𝑡
𝜕𝒓

f(l,m,n,i,j)=f(l,m,n,i,j) &-dfx(l,m,n,i,j)+dfx(l,m,n,i-1,j) &-dfy(l,m,n,i,j)+dfy(l,m,n,i,j-1)

流束の配列が大きい
線形移流
同じ計算方法 ⇒転置 ftmp(i,j,l,m,n)ftmp(i,j,l,m,n)=f(i,j,l,m,n) &“velocity_e” subroutine-dfx(i,j)+dfx(i-1,j) &-dfy(i,j)+dfy(i,j-1)
𝜕𝑓௦ 𝑞௦
𝜕𝑓௦
f(l,m,n,i,j)=f(l,m,n,i,j) &
+
𝑬ȉ
=0
𝜕𝑡 𝑚௦
𝜕𝒗
-dfx(l,m,n)+dfx(l-1,m,n) &-dfy(l,m,n)+dfy(l,m-1,n) &-dfz(l,m,n)+dfz(l,m,n-1)“velocity_b” subroutine

𝜕𝑓௦ 𝑞௦
𝜕𝑓௦
+
𝒗×𝑩 ȉ
=0
𝜕𝑡 𝑚௦
𝜕𝒗
回転
演算を各方向で分離(directional splitting)f(l,m,n,i,j)=f(l,m,n,i,j)-dfz(l,m,n)+dfz(l,m,n-1)f(l,m,n,i,j)=f(l,m,n,i,j)-dfy(l,m,n)+dfy(l,m-1,n)f(l,m,n,i,j)=f(l,m,n,i,j)-dfx(l,m,n)+dfx(l-1,m,n)f(l,m,n,i,j)=f(l,m,n,i,j)-dfx(l,m,n)+dfx(l-1,m,n)f(l,m,n,i,j)=f(l,m,n,i,j)-dfy(l,m,n)+dfy(l,m-1,n)f(l,m,n,i,j)=f(l,m,n,i,j)-dfz(l,m,n)+dfz(l,m,n-1)
6

プログラム２(as-is)性能分析● プロファイラ ループコスト*************************************************************************************
Application - loops*************************************************************************************
Cost
%
Nest
Kind
Exec
Start
End
----------------------------------------------------------------------------------------------89927
100.0000-----Application----------------------------------------------------------------------------------------------19623
21.82105DOOpenMP
137
157
velocity.velocity̲e̲5d.̲OMP̲1̲
137,221：flux累積和
4733
5.2632
5DO
OpenMP
251
258
velocity.velocity̲b̲5d.̲OMP̲2̲
4214
4.68605DOOpenMP
221
229
position.position̲5d.̲OMP̲2̲
3967
4.4114
5DO
OpenMP
126
133
velocity.velocity̲e̲5d.̲OMP̲1̲ 104,115,126,3525
3.91985DO
OpenMP
432
439
velocity.velocity̲b̲5d.̲OMP̲2̲ 251,286,396,432：
3428
3.81205DOOpenMP
157
159
position.position̲5d.̲OMP̲1̲
3350
3.7252
5DO
OpenMP
115
122
velocity.velocity̲e̲5d.̲OMP̲1̲ 関数pic5呼出3306
3.6763
5DO
OpenMP
286
293
velocity.velocity̲b̲5d.̲OMP̲2̲
3107
3.45505DOOpenMP
396
403
velocity.velocity̲b̲5d.̲OMP̲2̲ 157：転置3094
3.4406
5DO
OpenMP
104
111
velocity.velocity̲e̲5d.̲OMP̲1̲

velocity̲eは1タイム
ステップあたり2回呼ばれる

velocity̲e
137
5
138
5
139
5
140
5
141
5
142
5
143
5
144
5
145
5
146
5
147
5
148
5
149
5
150
5
151
5
152
5
153
5
154
5
155
5
156
5
157
5

*ループ長：42
oclによる手動unrolling
do ll=nvxs-1,nvxe+1
gfxy = abs((gfx(ll)*ffi(ll))* gfy(ll))*0.5d0gfyz = abs((gfy(ll)*ffi(ll))* gfz(ll))*0.5d0gfzx = abs((gfz(ll)*ffi(ll))* gfx(ll))*0.5d0gfxyz = abs((gfx(ll)*ffi(ll))*(gfy(ll)*ffi(ll))*gfz(ll))*inv3

p
p
p
p
p

10s
10v
10v
10v
10v

p
p
p
p

10v
10v
10v
10v

dfx(ll+l0,mm,nn) = dfx(ll+l0,mm,nn) +(gfx(ll)+(gfxyz-gfxy-gfzx)*sx)
dfx(ll+l0,mm,nn+nv) = dfx(ll+l0,mm,nn+nv)-(gfxyz-gfzx)*sx
dfx(ll+l0,mm+mv,nn) = dfx(ll+l0,mm+mv,nn)-(gfxyz-gfxy)*sx
dfx(ll+l0,mm+mv,nn+nv) = dfx(ll+l0,mm+mv,nn+nv)+ gfxyz
*sx

p
p
p
p

10m
10m
10m
10m

dfy(ll,mm+m0,nn) = dfy(ll,mm+m0,nn) +(gfy(ll)+(gfxyz-gfyz-gfxy)*sy)
dfy(ll+lv,mm+m0,nn) = dfy(ll+lv,mm+m0,nn)-(gfxyz-gfxy)*sydfy(ll,mm+m0,nn+nv) = dfy(ll,mm+m0,nn+nv)-(gfxyz-gfyz)*sydfy(ll+lv,mm+m0,nn+nv) = dfy(ll+lv,mm+m0,nn+nv)+ gfxyz
*sy

p
p
p
p
p

10m
10m
10m
10m
10v

dfz(ll,mm,nn+n0) = dfz(ll,mm,nn+n0) +(gfz(ll)+(gfxyz-gfzx-gfyz)*sz)
dfz(ll,mm+mv,nn+n0) = dfz(ll,mm+mv,nn+n0)-(gfxyz-gfyz)*sz
dfz(ll+lv,mm,nn+n0) = dfz(ll+lv,mm,nn+n0)-(gfxyz-gfzx)*sz
dfz(ll+lv,mm+mv,nn+n0) = dfz(ll+lv,mm+mv,nn+n0)+ gfxyz*sz
end doコンパイラによる
最適化なし

最適化阻害
7

プログラム３(ループ分割)性能分析
velocity̲e

129
130
131
132
133
134
135

5
5
5
5
5
5
5

p
p
p
p
p
p
p

138
139
140
141
142

5
5
5
5
5

p
p
p
p
p

145
146
147
148
149

5
5
5
5
5

p
p
p
p
p

152
153
154
155
156

159
160
161
162
163

5
5
5
5
5

5
5
5
5
5

p
p
p
p
p

p
p
p
p
p

<<<
<<<
<<<
<<<
10v
10v
10v
10v
10v
10v
10v
<<<
<<<
<<<
<<<
<<<
10v
10v
10v
10v
10v
<<<
<<<
10v
10v
10v
10v
10v
<<<
<<<
10v
10v
10v
10v
10v
<<<
<<<
10v
10v
10v
10v
10v

*ループ長：42 oclによる手動unrolling
SIMD(VL: 4)
SOFTWARE PIPELINING
PREFETCH(HARD) Expected by compiler :ff, (unknown), (unknown)
do ll=nvxs-1,nvxe+1
ffi = 1.0d0/ff(ll,mm,nn,ii,jj)
gfxy(ll) = abs((gfx(ll)*ffi)* gfy(ll))*0.5d0gfyz(ll) = abs((gfy(ll)*ffi)* gfz(ll))*0.5d0gfzx(ll) = abs((gfz(ll)*ffi)* gfx(ll))*0.5d0gfxyz(ll) = abs((gfx(ll)*ffi)*(gfy(ll)*ffi)*gfz(ll))*inv3
end doSIMD(VL: 4)
SOFTWARE PIPELINING
PREFETCH(SOFT) : 24
SEQUENTIAL : 24
(unknown): 24
do ll=nvxs-1,nvxe+1
dfx(ll+l0,mm,nn) = dfx(ll+l0,mm,nn) +(gfx(ll)+(gfxyz(ll)-gfxy(ll)-gfzx(ll))*sx)
dfy(ll,mm+m0,nn) = dfy(ll,mm+m0,nn) +(gfy(ll)+(gfxyz(ll)-gfyz(ll)-gfxy(ll))*sy)
dfz(ll,mm,nn+n0) = dfz(ll,mm,nn+n0) +(gfz(ll)+(gfxyz(ll)-gfzx(ll)-gfyz(ll))*sz)
end doSIMD(VL: 4)
SOFTWARE PIPELINING
do ll=nvxs-1,nvxe+1
dfx(ll+l0,mm,nn+nv) = dfx(ll+l0,mm,nn+nv)-(gfxyz(ll)-gfzx(ll))*sxdfy(ll+lv,mm+m0,nn) = dfy(ll+lv,mm+m0,nn)-(gfxyz(ll)-gfxy(ll))*sydfz(ll,mm+mv,nn+n0) = dfz(ll,mm+mv,nn+n0)-(gfxyz(ll)-gfyz(ll))*sz
end doSIMD(VL: 4)
SOFTWARE PIPELINING
do ll=nvxs-1,nvxe+1
dfx(ll+l0,mm+mv,nn) = dfx(ll+l0,mm+mv,nn)-(gfxyz(ll)-gfxy(ll))*sxdfy(ll,mm+m0,nn+nv) = dfy(ll,mm+m0,nn+nv)-(gfxyz(ll)-gfyz(ll))*sydfz(ll+lv,mm,nn+n0) = dfz(ll+lv,mm,nn+n0)-(gfxyz(ll)-gfzx(ll))*sz
end doSIMD(VL: 4)
SOFTWARE PIPELINING
do ll=nvxs-1,nvxe+1
dfx(ll+l0,mm+mv,nn+nv) = dfx(ll+l0,mm+mv,nn+nv)+ gfxyz(ll)*sxdfy(ll+lv,mm+m0,nn+nv) = dfy(ll+lv,mm+m0,nn+nv)+ gfxyz(ll)*sydfz(ll+lv,mm+mv,nn+n0) = dfz(ll+lv,mm+mv,nn+n0)+ gfxyz(ll)*sz
end doコンパイラによる最適化※loop A

ループ間のデータ受け渡しのために新たな1次元配列を使用
コンパイラによる最適化※loop B

コンパイラによる最適化※loop C

コンパイラによる最適化※loop D

コンパイラによる最適化※loop E

8プログラム３(ループ分割)性能分析● プロファイラ ループコスト：プルグラム３

velocity̲eは1タイム
ステップあたり2回呼ばれる*************************************************************************************
Application - loops*************************************************************************************
Cost
%
Nest
Kind
Exec
Start
End
----------------------------------------------------------------------------------------------77742
100.0000-----Application----------------------------------------------------------------------------------------------4730
6.0842
5DO
OpenMP
258
265
velocity.velocity̲b̲5d.̲OMP̲2̲
3601
4.63205DOOpenMP
119
126
velocity.velocity̲e̲5d.̲OMP̲1̲
3582
4.60755DO
OpenMP
439
446
velocity.velocity̲b̲5d.̲OMP̲2̲
3432
4.4146
5DO
OpenMP
293
300
velocity.velocity̲b̲5d.̲OMP̲2̲
3391
4.36195DO
OpenMP
157
159
position.position̲5d.̲OMP̲1̲
3356
4.31685DO
OpenMP
97
104
velocity.velocity̲e̲5d.̲OMP̲1̲
3237
4.16385DO
OpenMP
130
136
velocity.velocity̲e̲5d.̲OMP̲1̲
3227
4.15095DO
OpenMP
403
410
velocity.velocity̲b̲5d.̲OMP̲2̲
3091
3.97605DOOpenMP108
115
velocity.velocity̲e̲5d.̲OMP̲1̲
2849
3.6647
5DO
OpenMP
139
143
velocity.velocity̲e̲5d.̲OMP̲1̲

プログラム
position
velocity_e

As-is
(2)
1.29
3.19

ループ分割
(3)
1.03
2.48

97,108,119,
258,293,439,403：
関数pic5呼出157：転置130：loop A
139：loop B

比率
0.80
0.78

mpi_wtimeで計測した1タイムステップあたりの平均経過時間
velocity_bは変更なし

9分析区間velocity_e：プログラムの比較
 PAの実行時間とPA情報で比較
velocity_e
FX100(PA実行時間)

プログラム

スレッド数

As-is (2)ループ
分割 (3)

16

実行時間
(sec)

比率

15.68

1.00

12.34

0.79

 最も負荷が高いロープは、配列アクセスの依存関係で最適化がなされなかった。
 配列アクセスの依存性を無くすループ分割により最適化が促進され、約20%の計算時間の短縮に成功
 SIMD以外の最適化がかかっているループ(loopA, B)のコストが少し目立つ。
 Loop B-Eは配列インデックスの計算がある⇒!OCL NOSWPを指定したほうがよい？？

As-is
(2)

ループ
分割
(3)

0.79
4命令コミット

浮動小数点
演算待ち

浮動小数点ロード
メモリアクセス待ち

10

128reg/32regコンパイラ比較、FX100の測定条件
測定条件
項目

条件

コンパイラ

FX100用128regコンパイラVer.2.0.0 /ポストFX100向けチューニング用32regコンパイラ

翻訳時
オプション

-Kfast,ocl,simd=2 -x250

パラメータ

mpi_x = 1, mpi_y = 1, ns =2
nx = 128, ny = 64
nvx = 40, nvy = 40, nvz = 40
(実行サイズ約28GB)

スレッド数

16スレッド

実行条件

1ノード – ２プロセス

プログラム

ループ分割版(プログラム3)11

主要３区間の実行時間
 評価対象区間は、velocity_e、velocity_b、position の計3つのサブルーチン
レジスタ数
スレッド数
position
velocity_evelocity_b
計

128reg
16t
(秒)
5.40
12.39
11.67
29.46

32reg
16t
(秒)
6.27
17.97※
18.98
43.02

比率
1.16
1.45
1.62
1.46

velocity̲eは1タイムステッ
プあたり2回呼ばれるPAで計測した５タイムステップの経過時間※関数pic5呼出ループにおいて、!OCL UNROLLを外した場合（外さない場合は86.31秒：約7倍）
12

主要ループコスト● プロファイラ ループコスト 128reg16t（velocity̲eに!OCL UNROLLあり)*************************************************************************************
Application - loops*************************************************************************************
Cost
%
Nest
Kind
Exec
Start
End
----------------------------------------------------------------------------------------------90857
100.0000-----Application----------------------------------------------------------------------------------------------6613
7.27855DO
OpenMP
115
122
velocity.velocity̲e̲5d.̲OMP̲1̲
6212
6.8371
5DO
OpenMP
104
111
velocity.velocity̲e̲5d.̲OMP̲1̲
6076
6.6874
5DO
OpenMP
93
100
velocity.velocity̲e̲5d.̲OMP̲1̲
4341
4.77785DO
OpenMP
252
259
velocity.velocity̲b̲5d.̲OMP̲2̲
3897
4.2892
5DO
OpenMP
433
440
velocity.velocity̲b̲5d.̲OMP̲2̲
3642
4.00855DO
OpenMP
134
138
velocity.velocity̲e̲5d.̲OMP̲1̲
3423
3.76755DO
OpenMP
160
162
position.position̲5d.̲OMP̲1̲
3269
3.59805DOOpenMP
287
294
velocity.velocity̲b̲5d.̲OMP̲2̲
3145
3.46155DO
OpenMP
148
152
velocity.velocity̲e̲5d.̲OMP̲1̲
3091
3.40205DOOpenMP
141
145
velocity.velocity̲e̲5d.̲OMP̲1̲

● プロファイラ ループコスト

32reg16t（velocity̲eに!OCL UNROLLなし)*************************************************************************************
Application - loops*************************************************************************************
Cost
%
Nest
Kind
Exec
Start
End
----------------------------------------------------------------------------------------------115508
100.0000-----Application----------------------------------------------------------------------------------------------8260
7.15105DOOpenMP
115
122
velocity.velocity̲e̲5d.̲OMP̲1̲
7297
6.3173
5DO
OpenMP
93
100
velocity.velocity̲e̲5d.̲OMP̲1̲
7249
6.27585DO
OpenMP
104
111
velocity.velocity̲e̲5d.̲OMP̲1̲
5921
5.1261
5DO
OpenMP
252
259
velocity.velocity̲b̲5d.̲OMP̲2̲
4507
3.90195DO
OpenMP
397
404
velocity.velocity̲b̲5d.̲OMP̲2̲
4460
3.8612
5DO
OpenMP
360
367
velocity.velocity̲b̲5d.̲OMP̲2̲
4413
3.82055DO
OpenMP
287
294
velocity.velocity̲b̲5d.̲OMP̲2̲
4328
3.74695DO
OpenMP
433
440
velocity.velocity̲b̲5d.̲OMP̲2̲
4128
3.57385DO
OpenMP
323
330
velocity.velocity̲b̲5d.̲OMP̲2̲
3118
2.6994
5DO
OpenMP
160
162
position.position̲5d.̲OMP̲1̲

velocity̲eは
1タイムス
テップあたり2回呼ばれる*loop3*loop2*loop1*loop4*loop9*loop5loop4とloop9
loop5とloop8
loop6とloop7
は同じ構造*loop3*loop1*loop2*loop4*loop8*loop7*loop5*loop9*loop613

loop4とloop9
loop5とloop8
loop6とloop7
は同じ構造

高コストループ
velocity̲evelocity̲b

89
90

3
4

p
p

93
94
95
96
97
98
99
100

5
5
5
5
5
5
5
5

p
p
p
p
p
p
pi
p

v
v
v
v
v
v
v
v

104
105
106
107
108
109
110
111

5
5
5
5
5
5
5
5

p
p
p
p
p
p
pi
p

v
v
v
v
v
v
v
v

115
116
117
118
119
120
121
122

5
5
5
5
5
5
5
5

p
p
p
p
p
p
pi
p

v
v
v
v
v
v
v
v

do nn=nvzs-1,nvze+1
do mm=nvys-1,nvye+1*loop1 ⇒ ループ長：42
do ll=nvxs-1,nvxe+1
hm2=ff(ll-lv-lv,mm,nn,ii,jj)
hm1=ff(ll-lv,mm,nn,ii,jj)
hp0=ff(ll,mm,nn,ii,jj)
hp1=ff(ll+lv,mm,nn,ii,jj)
hp2=ff(ll+lv+lv,mm,nn,ii,jj)
gfx(ll)=pic5(hp2,hp1,hp0,hm1,hm2,fex)
end do*loop2 ⇒ ループ長：42
do ll=nvxs-1,nvxe+1
hm2=ff(ll,mm-mv-mv,nn,ii,jj)
hm1=ff(ll,mm-mv ,nn,ii,jj)
hp0=ff(ll,mm,nn,ii,jj)
hp1=ff(ll,mm+mv ,nn,ii,jj)
hp2=ff(ll,mm+mv+mv,nn,ii,jj)
gfy(ll)=pic5(hp2,hp1,hp0,hm1,hm2,fey)end do*loop3 ⇒ ループ長：42
do ll=nvxs-1,nvxe+1
hm2=ff(ll,mm,nn-nv-nv,ii,jj)
hm1=ff(ll,mm,nn-nv,ii,jj)hp0=ff(ll,mm,nn,ii,jj)
hp1=ff(ll,mm,nn+nv,ii,jj)
hp2=ff(ll,mm,nn+nv+nv,ii,jj)
gfz(ll)=pic5(hp2,hp1,hp0,hm1,hm2,fez)
end do関数pic5のインライン展開

関数pic5を呼び出すループが高コスト
形状はほぼ同じ。ループ長やメモリへのアクセスが異なる。

241
242

3
4

p
p

252
253
254
255
256
257
258
259
260

5
5
5
5
5
5
5
5
4

p 10v
p 10v
p 10v
p 10v
p 10v
p 10v
pi 10v
p 10v
p

276
277

3
4

p
p

287
288
289
290
291
292
293
294
295

5
5
5
5
5
5
5
5
4

p 10v
p 10v
p 10v
p 10v
p 10v
p 10v
pi 10v
p 10v
p

310
311

3
4

p
p

323
324
325
326
327
328
329
330
331

5
5
5
5
5
5
5
5
4

p 10v
p 10v
p 10v
p 10v
p 10v
p 10v
pi 10v
p 10v
p

do mm=nvys,nvye
do nn=nvzs-1,nvze*loop4 ⇒ ループ長：40 oclによる手動unrolling
do ll=nvxs,nvxe
hm2=ff(ll,mm,nn0(ll)-nnz(ll)-nnz(ll),ii,jj)
hm1=ff(ll,mm,nn0(ll)-nnz(ll),ii,jj)hp0=ff(ll,mm,nn0(ll),ii,jj)
hp1=ff(ll,mm,nn0(ll)+nnz(ll),ii,jj)
hp2=ff(ll,mm,nn0(ll)+nnz(ll)+nnz(ll),ii,jj)
dfz1(ll,nn)=pic5(hp2,hp1,hp0,hm1,hm2,vvz(ll))
end doend dodo nn=nvzs,nvze
do mm=nvys-1,nvye*loop5 ⇒ ループ長：40 oclによる手動unrolling
do ll=nvxs,nvxe
hm2=ff(ll,mm0(ll)-mmy(ll)-mmy(ll),nn,ii,jj)
hm1=ff(ll,mm0(ll)-mmy(ll),nn,ii,jj)
hp0=ff(ll,mm0(ll),nn,ii,jj)
hp1=ff(ll,mm0(ll)+mmy(ll),nn,ii,jj)
hp2=ff(ll,mm0(ll)+mmy(ll)+mmy(ll),nn,ii,jj)
dfy1(ll,mm)=pic5(hp2,hp1,hp0,hm1,hm2,vvy(ll))
end doend dodo nn=nvzs,nvze
do mm=nvys,nvye*loop6 ⇒ ループ長：40 oclによる手動unrolling
do ll=nvxs-1,nvxe
hm2=ff(ll0(ll)-llx(ll)-llx(ll),mm,nn,ii,jj)
hm1=ff(ll0(ll)-llx(ll),mm,nn,ii,jj)
hp0=ff(ll0(ll),mm,nn,ii,jj)
hp1=ff(ll0(ll)+llx(ll),mm,nn,ii,jj)
hp2=ff(ll0(ll)+llx(ll)+llx(ll),mm,nn,ii,jj)
dfx1(ll,mm)=pic5(hp2,hp1,hp0,hm1,hm2,vvx(ll))
end doend do関数pic5のインライン展開

14

主要行コスト● プロファイラ 行コスト 128reg16t*************************************************************************************
Application ‒ lines*************************************************************************************
Cost
%
Line
----------------------------------90857
100.0000-Application----------------------------------5491
6.0436121
velocity.velocity̲e̲5d.̲OMP̲1̲
5311
5.845499
velocity.velocity̲e̲5d.̲OMP̲1̲
5266
5.7959110velocity.velocity̲e̲5d.̲OMP̲1̲
2471
2.7197477
velocity.velocity̲b̲5d.̲OMP̲2̲
2443
2.6888477
velocity.velocity̲e̲5d.̲OMP̲1̲
2275
2.5039160
position.position̲5d.̲OMP̲1̲
2260
2.4874475
velocity.velocity̲b̲5d.̲OMP̲2̲
2091
2.3014276
position.position̲5d.̲OMP̲3̲
1916
2.1088166
velocity.velocity̲e̲5d.̲OMP̲1̲
1420
1.5629475
velocity.velocity̲e̲5d.̲OMP̲1̲

99,110,121:関数pic5呼出(1タイムステップあたり2回)● プロファイラ 行コスト

32reg16t*************************************************************************************
Application ‒ lines*************************************************************************************
Cost
%
Line
----------------------------------115508
100.0000-Application----------------------------------6605
5.7182121
velocity.velocity̲e̲5d.̲OMP̲1̲
6269
5.4273110velocity.velocity̲e̲5d.̲OMP̲1̲
6230
5.393699
velocity.velocity̲e̲5d.̲OMP̲1̲
4257
3.6855477
velocity.velocity̲b̲5d.̲OMP̲2̲
3280
2.8396475
velocity.velocity̲b̲5d.̲OMP̲2̲
2890
2.5020329
velocity.velocity̲b̲5d.̲OMP̲2̲
2830
2.4500366
velocity.velocity̲b̲5d.̲OMP̲2̲
2725
2.3591403
velocity.velocity̲b̲5d.̲OMP̲2̲
2692
2.3306477
velocity.velocity̲e̲5d.̲OMP̲1̲
2632
2.2786293
velocity.velocity̲b̲5d.̲OMP̲2̲

477,475：関数pic5内（velocity_bでは1タイプ
ステップあたり6回、velocity_eでは1タイプス
テップあたり3×2回）

293,329,366,403:関数pic5呼出(loop5,6,7,8)15

高コスト行(関数pic5内)関数pic5内
470
471
472
473
474
475
476
477

hmax1=max(max(hm1,hp0),min(hm1*2.0d0-hm2,hp0*2.0d0-hp1))
hmin1=min(min(hm1,hp0),max(hm1*2.0d0-hm2,hp0*2.0d0-hp1))
hmax2=max(max(hp1,hp0),min(hp0*2.0d0-hm1,hp1*2.0d0-hp2))
hmin2=min(min(hp1,hp0),max(hp0*2.0d0-hm1,hp1*2.0d0-hp2))
hmax=max(hmax1,hmax2)
hmin=max(gmin,min(hmin1,hmin2))

問題点、疑問点など
loop4とloop9のループコストが目立つが、128reg,32reg
コンパイラで共に、行コストには表れない。
loop4とloop9は同じ形状であるにもかかわらず、loop4のコストが高い。

関数pic5内で、470-473行目のほうが計算が重そうにもかかわらず、477,475行目のコストが128reg,32regコンパイラで共に目立つ。
32regコンパイラにおいて、477行目の性能劣化が顕著。
⇒ループ内で関数を呼び出さずに、関数内にループ計算する
ように構成変更すれば、コンパ
イラの最適化が促進？
16分析区間velocity_b：128reg と 32regコンパイラ の比較
 PAの実行時間とPA情報で比較
 128reg と 32regコンパイラの性能差は約57%で、レジスタ数減少の影響は大きい。
velocity_b

レジスタ数 スレッド数

FX100(PA実行時間)

128reg
32reg

16

実行時間
(sec)

比
率

11.67

1.00

18.98

1.62

FX100
128reg
16t

FX100
32reg
16t

1.62
4命令コミット レジスタ数減少の影響で、ループはアンローリングされているが、レジスタが少ないため、スピルが発生し、演算レイテンシ（浮動小数点演算待ち）が増加。
⇒ループアンローリング数の調整による性能改善の
可能性？？

浮動小数点
演算待ち

浮動小数点ロード
L1Dアクセス待ち

浮動小数点ロード
メモリアクセス待ち

17

ポストFX100、主要３区間の実行時間推定
 評価対象区間は、velocity_e、velocity_b、position の計3つのサブルーチン

レジスタ数
スレッド数
position
velocity_evelocity_b
計

FX100
PA実行
(a)
128reg
16t
(秒)
5.40
12.39
11.67
29.46

FX100
PA実行
(b)
32reg
16t
(秒)
6.27

ポスト
FX100
推定(c)
32reg
12t
(秒)
7.43

17.97※
18.98
43.02

23.31※
24.42
55.16

比率
c/b
1.185
1.297
1.287
1.282

 ポストFX100推定値は、32reg実行結果に対して約1.28倍 ⇒ 24スレッドなら約0.64倍
 演算ネックの区間は、実際の性能より低めの推定値が出る傾向にある※関数pic5呼出ループにおいて、!OCL UNROLLを外した場合 18

32regチュ−ニングのイメージ
As-is関数のベクトル化＋ループ分割

do nn=nvzs-1,nvze+1
do mm=nvys-1,nvye+1
do ll=nvxs-1,nvxe+1
hm2=ff(ll-lv-lv,mm,nn,ii,jj)
hm1=ff(ll-lv,mm,nn,ii,jj)
hp0=ff(ll,mm,nn,ii,jj)
hp1=ff(ll+lv,mm,nn,ii,jj)
hp2=ff(ll+lv+lv,mm,nn,ii,jj)
gfx(ll)=pic5(hp2,hp1,hp0,hm1,hm2,fex)
end do:::function pic5(hp2,hp1,hp0,hm1,hm2,vv)
real(kind=8), intent(in) :: hp2,hp1,hp0,hm1,hm2,vv
real(kind=8), intent(out):: pic5

do nn=nvzs-1,nvze+1
do mm=nvys-1,nvye+1
do ll=nvxs-1,nvxe+1
hm2(ll)=ff(ll-lv-lv,mm,nn,ii,jj)
hm1(ll)=ff(ll-lv,mm,nn,ii,jj)
hp0(ll)=ff(ll,mm,nn,ii,jj)
hp1(ll)=ff(ll+lv,mm,nn,ii,jj)
hp2(ll)=ff(ll+lv+lv,mm,nn,ii,jj)
ffx(ll)=fex
end docall pic(hp2,hp1,hp0,hm1,hm2,ffx,nvxe-nvxs+3,gfx):::
subroutine pic5(hp2,hp1,hp0,hm1,hm2,vv,nv,gf)
integer(kind=4) :: intent(in) :: nvreal(kind=8), intent(in) :: hp2(nv),hp1(nv),hp0(nv),hm1(nv),hm2(nv),vv(nv)real(kind=8), intent(out):: gf(nv)::pic5 =
end function

do ll=1,nv:
end dodo ll=1,nv

ループ分割により、velocity_bの実行時間の推定値が14.06秒、性能比c/bが0.74・c/aが1.2となった実績あり:end dodo ll=1,nv
gf(ll) =end doend subroutine

ループ分割による性能
チューニングは今後の課題※ループ分割に伴う作業配列の増加により性能が低下する可能性もある
19

4.3 プラズマ流体解析コード GT5D の Post-FX100 性能推定
日本原子力研究開発機構 井戸村 泰宏、伊奈 拓也
富士通株式会社 内藤 俊也、三吉 郁夫
4.3.1 は じ め に
核 融 合 プ ラ ズ マ 流 体 解 析 コ ー ド GT5D（ スライド p2、以降、ページ番号のみ表記する）
は磁場閉じ込め型の核融合炉における炉心プラズマの閉じ込め性能を決定する
乱 流 輸 送 現 象 の 評 価 を 目 的 と し て 開 発 さ れ た 。核融合プラズマの第一原理モデルは 5次元位相空間（3 次元空間×2 次元速度）におけるプラズマ粒子分布の発展を記述する移流・拡散方程式（ボルツマン方程式）で与えられる（p3）。核融合プラズマ流体解析ではこの方
程式を 5 次元格子上で離散化して CFD スキームを適用するため、3 次元流体モデルに比べて〜1002 倍程度大きい自由度の流体計算となる。
「京」を用いてようやく ITER 規模の炉心プラ
ズマの数値実験を断熱的電子モデルのイオン系乱流に対して行うことが可能になったが、運動論的電子モデルを含む電子系乱流の数値実験ではイオンに比べて熱速度が数十倍となる電子の運動を計算する必要があり、計算コストがさらに数十倍になるため、ポスト京のようなエクサスケールマシンが必要となる。しかしながら、現行機の SPARC64XIfx（FX100）
から将来の A64fx（ポスト京、あるいは、Post-FX100）への移行では、以下に示すような技
術的なギャップが存在する。
・SPARC から ARM へのアーキテクチャ変更
・レジスタ数の 128 から 32 への削減
・Out-of-Order 資源の増強
・SIMD 幅の 256bit から 512bit への増加
・プロセッサ（CMG）あたり演算コア数の 32（16）から 48（12）への変更
・メモリ機構の HMC（240+240GB/s）から HCM2（1024GB/s）への変更
このうち、本性能推定ではユーザプログラムの実装方法への影響が大きいと思われる、
レジスタ数の削減、および、メモリ機構の変更に関して、FX100 上で用意された 32 レジス
タ版富士通コンパイラ（32reg コンパイラ）と Post-FX100 性能推定ツールを用いて GT5D の主要カーネルであるクリロフソルバの性能推定を実施した（p1）。

4.3.2 プ ロ グ ラ ム 概 要
GT5D は 炉 心 プ ラ ズ マ 全 体 を 計 算 領 域 と し て プ ラ ズ マ 乱 流 や プ ラ ズ マ 分 布 の 大
域的な非線形発展を追跡する。この特徴を活かして、これまで非局所的な乱流輸送現象の研究や乱流輸送の装置サイズスケーリングの研究が行われてきた。
計算手法としては 5 次元移流・拡散方程式に差分法を適用し、時間積分に半陰解 法 を 用 い て い る （ p3） 。 移 流 項 を 陰 的 に 計 算 す る 差 分 ・ 陰 解 法 カ ー ネ ル は 共役 残 差 （ GCR） 法 、 あ る い は 、 省 通 信 一 般 化 最 小 残 差 （ CA-GMRES） 法 に よ っ て 実
装 さ れ て お り (p5)、コ ー ド の 主 要 コ ス ト を 占 め る（ p4）。並 列 化 手 法 と し て は M
PI と OpenMP を 用 い る ハ イ ブ リ ッ ド 並 列 モ デ ル を 採 用 し て い る 。 さ ら に 、 OpenM

P によって実装した通信スレッドによる通信隠蔽処理や階層的な多次元領域分
割 モ デ ル の 採 用 に よ り Oakforest-PACS 全 系 に 至 る 強 ス ケ ー リ ン グ を 達 成 し 、約1500 億 格 子 を 用 い る ITER 規 模 の 数 値 実 験 が 可 能 に な っ た （ p2） 。

4.3.3

128/32 レ ジ ス タ 版 コ ン パ イ ラ を 用 い た 性 能 測 定

p4 に 32reg/128reg コ ン パ イ ラ を 用 い て 測 定 し た GT5D の コ ス ト 分 布 を 示 す 。
こ こ で 、 問 題 サ イ ズ は (N x ,N y ,N z ,N v ,N w )=(320,320,32,96,20)=6x10 9 と し 、 MPI 並
列 数 は（ n x ,n y ,n w ）=(4,4,20)=320（ 160node x 2process x 16threads）と し た 。
GT5D は 移 流 項 の 4 次 元 差 分 計 算 、 衝 突 項 の 2 次 元 差 分 計 算 、 ポ ア ソ ン 方 程 式 の
2 次 元 有 限 要 素 法 、 1 次 元 FFT 等 、 異 な る 数 値 計 算 法 の カ ー ネ ル を 含 む が 、 多 くの 計 算 カ ー ネ ル で は 32reg コ ン パ イ ラ に よ る 性 能 劣 化 は 見 ら れ な か っ た 。一 方 、
全体コストの 8 割以上を占めるクリロフソルバ、および、移流項の差分計算カー ネ ル で は 大 き な 性 能 劣 化 が 見 ら れ た 。 こ こ で 、 ク リ ロ フ ソ ル バ の SpMV は 移 流
項 の 差 分 計 算 （ 非 対 称 行 列 、 17 点 ス テ ン シ ル ） で 与 え ら れ る こ と か ら 、 同 様 の性能劣化が発生しているものと推測される。以下では、さらにクリロフソルバ
内部のコスト分布を詳細に議論する。
p5 に GT5D で 実 装 さ れ て い る GCR 法 と CA-GMRES 法 の 演 算 特 性 を 示 す 。 GCR 法は SpMV カ ー ネ ル と 2 つ の DAXPY カ ー ネ ル で 構 成 さ れ 、 演 算 密 度 f/b=0.37 の 低
演 算 密 度 の ソ ル バ と な っ て い る 。一 方 、CA-GMRES 法 の 主 要 カ ー ネ ル は SpMV と 密
行 列 演 算 カ ー ネ ル （ DSYRK、 DGEMV） か ら 構 成 さ れ 、 f/b=1.38 の 高 演 算 密 度 の ソ
ル バ と な っ て い る 。 FX100 に お け る カ ー ネ ル 性 能 評 価 で は 、 DAXPY の 実 効 性 能 が
高 い GCR 法 の 実 効 性 能 が ル ー フ ラ イ ン モ デ ル の 85%と な る の に 対 し 、 CA-GMRES
法 の 実 効 性 能 は 61%と な る 。 し か し な が ら 、処 理 時 間 で 見 る と CA-GMRES 法 は GC
R 法 の 約 半 分 と な っ て い る 。 こ こ で 、 ル ー フ ラ イ ン モ デ ル は 各 環 境 に お け る str
eam triad 性 能 を 用 い た 。
p6 に GCR 法 、 CA-GMRES 法 の 各 カ ー ネ ル に 対 す る 32reg/128reg コ ン パ イ ラ 比較 を 示 す 。こ こ で 、問 題 サ イ ズ は (N x ,N y ,N z ,N v )=(160,160,32,96)と し 、MPI 並 列
数 は （ n x ,n y ） =(2,2)=4（ 2node x 2process x 16threads ） と し た 。 GCR 法 に お
い て は 、 ２ つ の DAXPY カ ー ネ ル の 32reg コ ン パ イ ラ に よ る 性 能 劣 化 は 見 ら れ な
い が 、SpMV カ ー ネ ル の 処 理 時 間 は 32reg コ ン パ イ ラ で 約 2 倍 に 増 大 す る 。一 方 、CA-GMRES 法 に お い て は 、 密 行 列 演 算 カ ー ネ ル の 性 能 劣 化 は 少 な い が 、 GCR 法 と共 通 の SpMV カ ー ネ ル に つ い て は 処 理 時 間 が 32reg コ ン パ イ ラ で 約 2 倍 に 増 大 す
る 。 こ こ で 、 GCR 法 と CA-GMRES 法 の SpMV カ ー ネ ル の 処 理 時 間 の 違 い は 、 GCR 法の SpMV カ ー ネ ル の み に 含 ま れ て い る ベ ク ト ル 内 積 演 算 処 理 に よ る も の で あ る 。
p7 に SpMV カ ー ネ ル の サ ン プ ル コ ー ド を 示 す 。こ の 差 分 計 算 は 4 次 元 移 流 項 の4 次精度中心差分（森西スキーム）で与えられるため、4 方向 5 点差分で合計 1
7 ス テ ン シ ル を 参 照 す る 。こ こ で 、ス テ ン シ ル 計 算 で 参 照 す る デ ー タ ア ク セ ス パ
タ ー ン は ス ト ラ イ ド ア ク セ ス と な っ て お り 、12MB の L2 キ ャ ッ シ ュ に 対 し て 最 内

3 重 ル ー プ の イ ン デ ッ ク ス 、 l、 k、 j ま で は オ ン キ ャ ッ シ ュ と な り 、 さ ら に 、 最
外 の i に つ い て も OpenMP に よ る dynamic ル ー プ 分 割 に よ っ て ス レ ッ ド 間 で キ ャ
ッ シ ュ 上 の デ ー タ 再 利 用 を 図 る こ と に よ り オ ン キ ャ ッ シ ュ と な る た め 、 L2 キ ャ
ッ シ ュ の 参 照 が 主 要 な コ ス ト を 占 め る [サイエンティフィックシステム研究会、ポスト
ペタアプリ性能評価 WG 成果報告書 2.8]。ま た 、こ の カ ー ネ ル は ル ー プ ボ デ ィ が 比 較
的大きく、レジスタ数の削減が処理性能に影響していると考えられる。この問題を回避するために、ユーザ側の最適化コードでは 4 次元差分の 4 重ループの最内の l 方向ループを 2 分割してループあたりのレジスタ数を削減した。ここ
で 、 分 割 点 は i,j 方 向 と k,l 方 向 、 各 8 ス テ ン シ ル ず つ に 分 割 し た 。 こ の ル ー
プ 分 割 に よ り 変 数 df へ の ロ ー ド /ス ト ア が 1 回 ず つ 増 え る が 、 最 内 ル ー プ の イ
ンデックス l は 4 次元配列の最内、すなわち、連続アクセスなので、増大する
ロ ー ド /ス ト ア は キ ャ ッ シ ュ 上 で 処 理 さ れ 、そ の ペ ナ ル テ ィ は 小 さ い と 考 え ら れる 。p6 に 最 適 化 後 の SpMV カ ー ネ ル の 処 理 性 能 を 示 す が 、処 理 時 間 の 増 大 は １ 割
以 下 に 抑 制 さ れ て お り 、 32reg/128reg コ ン パ イ ラ で ほ ぼ 同 等 の 性 能 が 得 ら れ る
こ と が わ か っ た 。こ の 修 正 カ ー ネ ル を 用 い て 32 レ ジ ス タ の Intel 環 境 で も 性 能
比 較 を 行 っ た が 、Intel 環 境 で は 最 適 化 前 後 で 性 能 に 違 い が 見 ら れ な か っ た 。これ は 、Intel コ ン パ イ ラ で は こ の よ う な 最 適 化 が 自 動 で 行 わ れ て い る こ と を 示 唆
している。

4.3.4

SpMV カ ー ネ ル の 詳 細 分 析 、 最 適 化

p9 で は 、 性 能 劣 化 が 見 ら れ た SpMV カ ー ネ ル に つ い て 詳 細 PA を 用 い た 分 析 を
実 施 し た 。 上 記 の 事 例 と 同 様 に 32reg コ ン パ イ ラ で 性 能 劣 化 が 見 ら れ る が 、 この 例 で は 、 後 で 議 論 す る Post-FX100 推 定 用 に 12 ス レ ッ ド で も 割 り 切 れ る よ うに 問 題 サ イ ズ を (N x ,N y ,N z ,N v )=(48,48,16,128)の 逐 次 処 理（ 16thread、も し く は 、
12thread） と し た こ と 、 お よ び 、 富 士 通 側 で コ ン パ イ ラ オ プ シ ョ ン を 調 整 し た
こ と の 影 響 に よ り 、 32reg/128reg 性 能 差 は 16%程 度 に 抑 制 さ れ て い る 。 PA 情 報
か ら 、 128reg コ ン パ イ ラ で は SpMV カ ー ネ ル の L2 ス ル ー プ ッ ト が 大 き い が 、 32
reg コ ン パ イ ラ で は L2 ス ル ー プ ッ ト が 約 15%低 減 し 、ロ ー ド /ス ト ア 命 令 数 が 約1 割増加している。これはレジスタスピルによる命令数の増加を示唆している。
こ の 問 題 の 回 避 に 向 け て 、ル ー プ 分 割 に よ る 最 適 化 を 実 施 し た 。p10、11 に p
7 に 示 す SpMV カ ー ネ ル の 最 内 2 重 ル ー プ を 切 り 出 し た サ ン プ ル コ ー ド を 示 す が 、
富 士 通 側 の 最 適 化 で は 、分 割 点 を i 方 向 と j,k,l 方 向 、そ れ ぞ れ 、4 ス テ ン シ ル 、
12 ス テ ン シ ル と 非 対 称 に 分 割 し て い る 。 ま た 、 ル ー プ 分 割 の 次 元 に 関 し て も 最
内 2 次 元 ル ー プ を 対 象 と し た 。 さ ら に 、 比 較 的 高 い L2 ミ ス dm 率 が 示 唆 す る メ
モ リ レ イ テ ン シ を 隠 蔽 す る た め に 、ル ー プ 分 割 後 の 各 2 次 元 ル ー プ に 対 し て 、unroll 数 を 2、 4 と 指 定 し て ソ フ ト ウ ェ ア プ リ フ ェ ッ チ を 手 動 で 指 定 し た 。 こ の結 果 、 p12 に 示 す よ う に 、 128reg コ ン パ イ ラ と ほ ど 同 等 の 処 理 時 間 が 得 ら れ 、 L
2 スループットも回復した。

4.3.5

GCR ソ ル バ の Post-FX100 性 能 推 定

最適化後の GCR ソルバを用いて Post-FX100 性能推定を実施し、2 つの DAXPY カーネルとSpMV カーネルの処理性能を比較した。ここでは FX100、Post-FX100 それぞれの 1CMG を用いた処理性能を比較した。p1 に各環境のスペックを示すが、FX100 の 1CMG は 16core、0.5TF
lops、メモリ帯域幅 120+120GB/s（in/out）、メモリ B/F=0.48、L2 キャッシュ B/F=2.0、L1
キャッシュ B/F=4.0、L2 キャッシュ 12MB となるのに対し、Post-FX100 の 1CMG は 12core、
0.675TFlops、256GB/s、メモリ B/F=0.38、L2 キャッシュ B/F=1.3、L1 キャッシュ B/F=4.0、
L2 キャッシュ 8MB となる。このため FX100/Post-FX100 理論性能比は演算性能で 1.35 倍、
メモリバンド幅で 1.06 倍となる。性能推定ツールでは FX100 の PA 情報に基づいて Post-F
X100 性能推定を行う。演算に関しては、SIMD 化されている命令については命令数を半分にカウントしている。メモリに関しては、メモリアーキテクチャの違いを反映してロード/ス
トアの処理コストや L1/L2 キャッシュアクセスのコストを換算しているがが、キャッシュ
サイズの違いは考慮していない。
p14、15 に DAXPY カーネルの性能を示す。DAXPY カーネルでは性能劣化が見られなかったことから性能評価は asis カーネルで実施した。2 つの DAXPY カーネルでほぼ同じ処理特性が得られており、どちらもメモリスループットがボトルネックとなっている。推定結果における Post-FX100 における 1CMG あたりの処理時間は FX100 の 0.75 倍となっている。一方、ルーフラインモデルで評価した理論性能比は DAXPY1（f/b=0.08）で 0.93 倍、DAXPY2（f/b=0.21）で 0.92 倍となる。ここで、ルーフラインモデルは上記理論性能を用いた。この推定結果は理論性能を上回る結果となっているが、これには Post-FX100 におけるメモリアー
キテクチャの変更が大きく影響している。FX100 の HMC がロード/ストアそれぞれ 1:1 の設計となっており、ロード 2:ストア 1 となる stream triad の実効メモリバンド幅が 6 割程度にとどまっていた[サイエンティフィックシステム研究会、ポストペタアプリ性能評価 WG
成果報告書 2.1]のに対し、Post-FX100 の HBM2 では非対称なロード/ストアでもメモリバン
ド幅をフルに活用でき、stream triad の実効メモリバンド幅が約 8 割に向上している[http://pr.fujitsu.com/jp/news/2018/08/22-1.html]。
p16 に SpMV カーネルの性能を示す。カーネルの演算特性としては、FX100 と同様に L2 ス
ループットがボトルネックとなっているが、SIMD 幅の拡大を反映して水色や青で示す命令
コミットの処理時間が削減されている。推定結果における Post-FX100 における 1CMG あたりの処理時間は FX100 とほぼ同じ処理時間となっている。ルーフラインモデルに基づく理論性能比では SpMV（f/b=1.03）に対して 0.87 倍となることから、この推定結果は理論性能を下回る結果となっている。このカーネルは L2 ビジー時間が主要なコストを占めることから、実行性能の劣化は L2 キャッシュ帯域幅が FX100 の B/F=2.0（1TB/s）から Post-FX100の B/F=1.3（877.5GB/s）に低下したことが影響していると考えられる。

4.3.6

まとめ

プラズマ流体解析コード GT5D を対象として Post-FX100 の性能推定を行った。GT5D コー
ドは移 流 項 の 4 次 元 差 分 計 算 、衝 突 項 の 2 次 元 差 分 計 算 、ポ ア ソ ン 方 程 式 の 2 次
元 有 限 要 素 法 、 1 次 元 FFT 等 、多彩な演算カーネルを含むコードであるが、4 次元差分計算以外は 32reg コンパイラによる性能劣化の影響は少ないことがわかった。これは、多くの既存アプリで大幅なコード修正なしに Post-FX100 の性能を引き出せる可能性を示唆している。実際に、GT5D 以外にも原子力機構の 3 次元多相多成分熱流動解析コードの評価を行ったところ、全てのカーネルで性能劣化が見られなかった。
4 次元差分計算の SpMV カーネルでは PA 情報からレジスタスピルによるものと推定される性能劣化が発生していたが、ループ分割による最適化を適用することで性能劣化を回避することに成功した。ただし、この実装に関してはユーザ側最適化と富士通側最適化でループ分割を行う次元や分割点に関して異なる実装が行われ、さらに、富士通側最適化では
分割後のループに対する unroll、prefech 等の手動最適化が適用された。これらの最適化
方法には任意性があるが、このような最適化の自動処理に向けた Post-FX100 コンパイラの開発が進展しており、今後、Post-FX100 の使い勝手が大幅に向上すると期待される。
Post-FX100 の性能推定では GCR ソルバの SpMV カーネルと 2 つの DAXPY カーネルの性能推定を実施した。DAXPY カーネルに関してはメモリアーキテクチャの変更により非対称なロード/ストア負荷のカーネルに対して性能向上が期待できることがわかった。一方、ほとんどのデータがオンキャッシュで処理され L2 スループットの負荷が大きい SpMV カーネルに関しては 1CMG あたりの処理時間が FX100 と Post-FX100 で同程度となり、ルーフラインモデ
ルの理論性能予測を下回る推定結果となった。この性能劣化は L2 キャッシュ帯域幅の低下によるものと考えられる。
今回の性能推定は FX100 の PA 情報に基づく性能推定ツールを用いたため、Out-of-Order
機構や ARM-SVE の新たな命令セットの効果等を見ることができなかった。これらの課題については、今後、ポスト京シミュレータを用いてさらに精度の高い性能推定を行うことが必要である。

性能評価環境

測定環境：名大FX100+128reg/32regコンパイラ
SPARC64XIfx: 1TFlops, 240+240GB/s, 16x2core, 256bitSIMD, HMC

評価単位CMG: 0.5TFlops, 240GB/s, 16corempifrtpx -X9 -Kfast -Kparallel -Karray_private -Kpreex -Kopenmp -Kocl -Kmfunc=2 Kprefetch_cache_level=all -CcdRR8 -Cpp -Koptmsg=2 -Y0,/center/fjlocal/opt/Simulation_compiler/bin -YI,/center/fjlocal/opt/Simulation_compiler/include
推定対象：ポストFX100
A64FX: 2.7TFlops, 256x4GB/s, 12x4core, 512bitSIMD, HBM2

評価単位CMG: 0.675TFlops, 256GB/s, 12core
1プラズマ流体解析コードGT5Dの概要







Gyrokinetic Toroidal 5D Eulerian code
5次元ジャイロ運動論モデルによる核融合炉心プラズマの第一原理乱流計算を行い、乱流輸送特性を評価
計算手法
MPI+OpenMPハイブリッド並列モデル
5次元位相空間の移流・拡散→差分法
時間積分→半陰的ルンゲ・クッタ法
差分・陰解法カーネル→省通信クリロフ部分空間法+通信隠蔽計算規模（ITER規模）
768x64x768x128x32 ~ 1500億格子~105ステップ→約5Mnode-h/ショット
参考文献
Idomura et al., Comput.Phys.Commun.2008
Idomura et al., Nucl. Fusion 2009
Idomura et al., Int. J. HPC Appl. 2014
Idomura, J. Comput. Phys. 2016
Idomura et al., ScalA17@SC17

5次元位相空間（速度2次元）の摂動粒子分布

v^
v//
v^

z
y

v//

x

Oakforest-PACSにおけるGT5D@ITER規模の強スケーリング
50000
40000

msec/step



Other
Field
Krylov(Calc)
Krylov(Halo)

Collision
Nonlinear
Krylov(Wait)
Krylov(Reduce)

30000
20000
10000
0
CA-GMRES/KNL
8x8x32=2048

CA-GMRES/KNL
8x16x32=4096

CA-GMRES/KNL
16x16x32=8192

計算モデル[Idomura,JCP07,CPC08,NF09,PFR14,JCP16]

プラズマ分布の移流・拡散方程式

熱運動や乱流場による移流

クーロン衝突効果

中心差分（4次精度森西）+半陰的時間積分（2次精度加法半陰的R-K）
線形移流項の差分・陰解法の非対称行列に対するクリロフソルバ

乱流場を決めるポアソン方程式

Z

z
a

y



2次元有限要素法
FFT

q
j

R

R0
3

GT5Dの全体コスト比較
問題サイズ：NxxNyxNzxNvxNw=320x320x32x96x20=6x109
並列化：nxxnyxnw=4x4x20=160 (160node x 2process x 16threads)
ms/step

128reg32reg

32reg/128reg

dn （密度）

20

22

1.09

field （電場）

247

247

1.00

drift （移流速度）

93

98

1.05

l4d_nl （非線形移流）

826

1131

1.37

l4d_l （線形移流）

466

806

1.73

rk （時間積分）

197

198

1.00

bc （境界条件）

922

1289

1.40

col （粒子衝突）

125

138

1.10

GCR （クリロフソルバ）

22966

29196

1.27

total

25862

33124

1.28

クリロフソルバが全体コストの8割以上を占める
クリロフソルバ、および、移流項の差分計算カーネルで顕著な性能劣化
他カーネルではあまり性能劣化は見られない
4

GT5Dにおけるクリロフソルバの演算特性[Idomura,ScalA17@SC17]

行列特性：非対称行列,17点ステンシル（5点差分✕4次元）
問題サイズ：NxxNyxNzxNv=160x160x32x96
並列化：nxxny=2x1 or 2x2 (2node x 2 or 1process)f/b:演算密度, t:msec/iteration, P:GFLOPS, RRL:Roofline比, RICEX:ICEX比

一般化共役残差法GCR

差分計算（SpMV）+ベクトル演算（AXPY）

低演算密度

FX100ではルーフライン理論性能の約85%の実効性能
省通信一般化最小残差法CA-GMRES

差分計算（SpMV）+密行列演算（SYRK）+行列ベクトル積（GEMV）

高演算密度

FX100ではルーフライン理論性能の約60%の実効性能

532reg/128regでのクリロフソルバ性能比較
GCR
msec/iteration

128reg32reg

32reg/128reg

32reg(opt)

32reg/128reg(opt)

SpMV1

7.90

16.38

2.07

8.68

1.10

DAXPY1

7.22

7.22

1.00

7.22

1.00

DAXPY2

7.07

7.06

1.00

7.05

1.00

TOTAL

22.19

30.66

1.38

22.95

1.03

msec/iteration

128reg32reg

32reg/128reg

32reg(opt)

32reg/128reg(opt)

SpMV2

6.55

13.21

2.02

6.80

1.04

DSYRK

3.37

3.95

1.17

3.95

1.17

DGEMV

1.93

1.89

0.98

1.92

1.00

TOTAL

11.84

19.05

1.61

12.67

1.07

CA-GMRESDAXPY、密行列演算（DSYRK,DGEMV）では性能劣化はあまり見られない

SpMV（差分計算）で大きい性能劣化が見られる
 ループ分割によるループボディ調整で性能を改善→32regによる性能劣化はOut of Order機能なしにほぼ解消？
6ユーザ側で実施した32reg向けSpMV最適化
オリジナル

最適化

real*8 f(-1:nv+2,-1:nz+2,-1:ny+2,-1:nx+2)!$OMP DO COLAPSE(2) SCHEDULE(dynamic,1)
do i = 1,nx (nx=80)
do j = 1,ny (ny=80)
compute vxyv~vxl2 from vx(l,j,i),vy(l,j,i),vz(l,j,i),vv(l,j,i)do k = 1,nz (nz=32)
do l = 1,nv (nv=96)
df(l,k,j,i)=&vxyv(l)*f(l,k,j,i)&+ vvr(l) *f(l+1,k,j,i) - vvl(l) *f(l-1,k,j,i)&- vvr2(l)*f(l+2,k,j,i) + vvl2(l)*f(l-2,k,j,i)&+vzl2(l)*(f(l,k-2,j,i)-f(l,k+2,j,i))&+vzl(l) *(f(l,k+1,j,i)-f(l,k-1,j,i))&+ vyr(l) *f(l,k,j+1,i) - vyl(l) *f(l,k,j-1,i)&- vyr2(l)*f(l,k,j+2,i) + vyl2(l)*f(l,k,j-2,i)&+ vxr(l) *f(l,k,j,i+1) - vxl(l) *f(l,k,j,i-1)&- vxr2(l)*f(l,k,j,i+2) + vxl2(l)*f(l,k,j,i-2)
enddo
enddo
enddo
enddo

real*8 f(-1:nv+2,-1:nz+2,-1:ny+2,-1:nx+2)!$OMP DO COLAPSE(2) SCHEDULE(dynamic,1)
do i = 1,nx (nx=80)
do j = 1,ny (ny=80)
compute vxyv~vxl2 from vx(l,j,i),vy(l,j,i),vz(l,j,i),vv(l,j,i)do k = 1,nz (nz=32)
do l = 1,nv (nv=96)
df(l,k,j,i)=&vxyv(l)*f(l,k,j,i)&+ vvr(l) *f(l+1,k,j,i) - vvl(l) *f(l-1,k,j,i)&- vvr2(l)*f(l+2,k,j,i) + vvl2(l)*f(l-2,k,j,i)&+vzl2(l)*(f(l,k-2,j,i)-f(l,k+2,j,i))&+vzl(l) *(f(l,k+1,j,i)-f(l,k-1,j,i))
enddo
do l = 1,nv (nv=96)
df(l,k,j,i) = df(l,k,j,i)+&+ vyr(l) *f(l,k,j+1,i) - vyl(l) *f(l,k,j-1,i)&- vyr2(l)*f(l,k,j+2,i) + vyl2(l)*f(l,k,j-2,i)&+ vxr(l) *f(l,k,j,i+1) - vxl(l) *f(l,k,j,i-1)&- vxr2(l)*f(l,k,j,i+2) + vxl2(l)*f(l,k,j,i-2)
enddo
enddo
enddo
enddo

この最適化はIntel環境（32reg）の性能には影響しない
7

SpMVカーネルの分析、最適化
 GCRカーネル
行列特性：非対称行列,17点ステンシル（5点差分✕4次元）
問題サイズ：NxxNyxNzxNv=48x48x16x128（16と12で割り切れるようにNxxNy=48x48と調整）
コンパイラオプション：(共通)-xpre̲diff -Kfast -Kopenmp -Kocl,optmsg=2 -V -X9 -Kmfunc=2 -Cpp -Kloop̲nofusion-CcdRR8 -Kswp̲strong -Kprefetch̲nostrong -Nlst=t,lst=i,lst=x -Kprefetch̲sequential=soft(DAXPY) -Kprefetch̲cache̲level=all(SpMV)-Kprefetch̲cache̲level=1

FUJITSU CONFIDENTIAL

8

Copyright 2018 FUJITSU LIMITED128reg と 32regコンパイラ の比較
 PAの実行時間で比較
 128reg と 32regコンパイラで、性能差が16%程度出ており、レジスタ数の影響は、やや大きい。
DIFF

コード

FX100 (PA実行時間)

レジスタ数

スレッド数

128regasis

実行時間(sec)

比率

0.37

1.00

0.43

1.16

16

32reg

 PA情報で比較
 128regはL2ビジー率が高めである。
 128reg,32reg共にL2デマンドミス率が高い。
 32regは、 128regに比べてスピル・フィル命令数が増加し、命令コミットが増加している。
 128reg、32reg共にSWPLが適用される。
DIFF

コード

レジスタ数

FX100
(PA情報)

asis

128reg

asis

32reg

DIFF

コード

レジスタ数

FX100
(PA情報)

asis

128reg

asis

32reg

FUJITSU CONFIDENTIAL

スレッド数

16
スレッド数

16

実行時間
(sec)

L1
ビジー率

L2
ビジー率

L2
スループット(GB/sec)

メモリ
ビジー率

メモリ
スループット(GB/sec)

0.37

69%

67%

263.09

69%

96.82

0.43

75%

57%

225.21

58%

80.64

ロード・
ストア数

L2
ミス率

L2 ミス
dm 率

2.79E+10 1.35% 16.83%

0.39%

44.34%

3.08E+10 1.22% 14.52%

0.35%

50.10%

9

L1D
ミス率

L1D ミスdm 率

Copyright 2018 FUJITSU LIMITEDチューニング：ループ分割位置
 ループ分割位置について (2次元-2分割)
2020 3 p
do k = 1,nz
⇒ nz は 162021 4 p 2v
do l = 1,nv
⇒ nv は 128
2023 4 p 2v
df(l,k,j,i)=2024 4&vxyv(l)*f(l,k,j,i)
2025 4&+ vvr(l) *f(l+1,k,j,i)2026 4&- vvl(l) *f(l-1,k,j,i)2027 4&- vvr2(l)*f(l+2,k,j,i)2028 4&+ vvl2(l)*f(l-2,k,j,i)2029 4&+ vzl2(l)*(f(l,k-2,j,i)-f(l,k+2,j,i))
2030 4&+ vzl(l) *(f(l,k+1,j,i)-f(l,k-1,j,i))
2031 4&+ vyr(l) *f(l,k,j+1,i)
2032 4&- vyl(l) *f(l,k,j-1,i)
2033 4&- vyr2(l)*f(l,k,j+2,i)2034 4&+ vyl2(l)*f(l,k,j-2,i)!--------------------------------------------------------------!
2035 4&+ vxr(l) *f(l,k,j,i+1)2036 4&- vxl(l) *f(l,k,j,i-1)
2037 4&- vxr2(l)*f(l,k,j,i+2)2038 4&+ vxl2(l)*f(l,k,j,i-2)
2040 4 p 2v
dfc̲(l)=dfc̲(l)+df(l,k,j,i)*fc(l,k,j,i)
2041 4 p 2v
enddo
2042 3 p
enddo

FUJITSU CONFIDENTIAL

2020
2021
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034

do k = 1,nz
⇒ nz は 16分
do l = 1,nv
⇒ nv は 128
割
df(l,k,j,i)=1&
vxyv(l)*f(l,k,j,i)&+ vvr(l) *f(l+1,k,j,i)&- vvl(l) *f(l-1,k,j,i)&- vvr2(l)*f(l+2,k,j,i)&+ vvl2(l)*f(l-2,k,j,i)&+ vzl2(l)*(f(l,k-2,j,i)-f(l,k+2,j,i))&+ vzl(l) *(f(l,k+1,j,i)-f(l,k-1,j,i))&+ vyr(l) *f(l,k,j+1,i)&- vyl(l) *f(l,k,j-1,i)&- vyr2(l)*f(l,k,j+2,i)&+ vyl2(l)*f(l,k,j-2,i)
enddo
enddo
!--------------------------------------------------------------!
do k = 1,nz
⇒ nz は 16分
do l = 1,nv
⇒ nv は 128
割2035 4&+ vxr(l) *f(l,k,j,i+1)22036 4&- vxl(l) *f(l,k,j,i-1)
2037 4&- vxr2(l)*f(l,k,j,i+2)2038 4&+ vxl2(l)*f(l,k,j,i-2)
2040 4 p 2v
dfc̲(l)=dfc̲(l)+df(l,k,j,i)*fc(l,k,j,i)
2041 4 p 2v
enddo
2042 3 p
enddo

10

3
4
4
4
4
4
4
4
4
4
4
4
4
4

p
p
p

2v
2v

Copyright 2018 FUJITSU LIMITEDチューニング：最適化の状況
 最適化の状況 (DIFF 区間) (ループ2次元-2分割 + 手動L2プリフェッチ・チューニング)

2040
2041
2042
2043
2044
2045
2047
2048
2049
2050
2051
2053
2067

3
3
3
3
3
3
3
3
3
3
3
3
3

p

2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083

4
4
4
4
4
4
4
4
4
4
4
4
4
4
3
2

p
p

2v
2v

p
p

2v

p
p
p
p
p
p
p
p

<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
PREFETCH: 8<<<
f: 4, df: 4<<< Loop-information End >>>do k = 1,nz
!ocl noprefetch!ocl prefetch̲read(f(1,k,j+3,i),level=2,strong=1)!ocl prefetch̲read(f(33,k,j+3,i),level=2,strong=1)!ocl prefetch̲read(f(65,k,j+3,i),level=2,strong=1)!ocl prefetch̲read(f(97,k,j+3,i),level=2,strong=1)!---------------------------------------------------!ocl prefetch̲write(df(1,k,j+1,i),level=2,strong=1)!ocl prefetch̲write(df(33,k,j+1,i),level=2,strong=1)!ocl prefetch̲write(df(65,k,j+1,i),level=2,strong=1)!ocl prefetch̲write(df(97,k,j+1,i),level=2,strong=1)!---------------------------------------------------!ocl unroll(2)<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)<<<
SOFTWARE PIPELINING<<< Loop-information End >>>do l = 1,nv
df(l,k,j,i)=&vxyv(l)*f(l,k,j,i)&+ vvr(l) *f(l+1,k,j,i)&- vvl(l) *f(l-1,k,j,i)&- vvr2(l)*f(l+2,k,j,i)&+ vvl2(l)*f(l-2,k,j,i)&+ vzl2(l)*(f(l,k-2,j,i)-f(l,k+2,j,i))&+ vzl(l) *(f(l,k+1,j,i)-f(l,k-1,j,i))&+ vyr(l) *f(l,k,j+1,i)&- vyl(l) *f(l,k,j-1,i)&- vyr2(l)*f(l,k,j+2,i)&+ vyl2(l)*f(l,k,j-2,i)
enddo
enddo
!-----------<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
PREFETCH: 24<<<
f: 16, df: 4, fc: 4<<< Loop-information End >>>FUJITSU CONFIDENTIAL

分
割
1

11

2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2158

3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3
3

p

2159
2160
2161
2162
2163
2164
2166
2167
2168

4
4
4
4
4
4
4
4
3

p
p

4v
4v

p
p
p

4v
4v

p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p
p

do k = 1,nz
!ocl noprefetch!ocl prefetch̲read(f(1,k,j+1,i-2),level=2,strong=1)!ocl prefetch̲read(f(1,k,j+1,i-1),level=2,strong=1)!ocl prefetch̲read(f(1,k,j+1,i+1),level=2,strong=1)!ocl prefetch̲read(f(1,k,j+1,i+2),level=2,strong=1)!ocl prefetch̲write(df(1,k,j+1,i),level=2,strong=1)!ocl prefetch̲read(fc(1,k,j+1,i),level=2,strong=1)!---------------------------------------------------!ocl prefetch̲read(f(33,k,j+1,i-2),level=2,strong=1)!ocl prefetch̲read(f(33,k,j+1,i-1),level=2,strong=1)!ocl prefetch̲read(f(33,k,j+1,i+1),level=2,strong=1)!ocl prefetch̲read(f(33,k,j+1,i+2),level=2,strong=1)!ocl prefetch̲write(df(33,k,j+1,i),level=2,strong=1)!ocl prefetch̲read(fc(33,k,j+1,i),level=2,strong=1)!---------------------------------------------------!ocl prefetch̲read(f(65,k,j+1,i-2),level=2,strong=1)!ocl prefetch̲read(f(65,k,j+1,i-1),level=2,strong=1)!ocl prefetch̲read(f(65,k,j+1,i+1),level=2,strong=1)!ocl prefetch̲read(f(65,k,j+1,i+2),level=2,strong=1)!ocl prefetch̲write(df(65,k,j+1,i),level=2,strong=1)!ocl prefetch̲read(fc(65,k,j+1,i),level=2,strong=1)!---------------------------------------------------!ocl prefetch̲read(f(97,k,j+1,i-2),level=2,strong=1)!ocl prefetch̲read(f(97,k,j+1,i-1),level=2,strong=1)!ocl prefetch̲read(f(97,k,j+1,i+1),level=2,strong=1)!ocl prefetch̲read(f(97,k,j+1,i+2),level=2,strong=1)!ocl prefetch̲write(df(97,k,j+1,i),level=2,strong=1)!ocl prefetch̲read(fc(97,k,j+1,i),level=2,strong=1)!---------------------------------------------------!ocl unroll(4)<<< Loop-information Start >>><<< [OPTIMIZATION]<<<
SIMD(VL: 4)<<<
SOFTWARE PIPELINING<<< Loop-information End >>>do l = 1,nv
df(l,k,j,i)= df(l,k,j,i)&+ vxr(l) *f(l,k,j,i+1)&- vxl(l) *f(l,k,j,i-1)&- vxr2(l)*f(l,k,j,i+2)&+ vxl2(l)*f(l,k,j,i-2)
dfc̲(l)=dfc̲(l)+df(l,k,j,i)*fc(l,k,j,i)
enddo
enddo

分
割
2

Copyright 2018 FUJITSU LIMITEDチューニングによる性能改善
 PA情報で比較 (チューニングの効果)
 ループ2分割 ＋ L2プリフェッチ によるチューニングの実施
 2次元で2分割した各ループの unroll 数を調整(2と4)して SWPL を促進 1ループ中の処理を小さくすることにより、スピル・フィルを削減
 L2デマンドミス率が高いため、メモリレイテンシの隠蔽で、L2プリフェッチを実施
DIFF
FX100
(PA情報)

レジスタ数

asis

128reg

asis
⇒tune

DIFF
FX100
(PA情報)

コード

32reg

コード

レジスタ数

asis

128reg

asis
⇒tune

32reg

スレッド数

16

スレッド数

16

実行時間
(sec)

L1
ビジー率

L2
ビジー率

L2
スループット(GB/sec)

メモリ
ビジー率

メモリ
スループット(GB/sec)

0.37

69%

67%

263.09

69%

96.82

0.43

75%

57%

225.21

58%

80.64

0.38

71%

68%

261.57

72%

99.65

ロード・
ストア数

L2
ミス率

L2 ミス
dm 率

2.79E+10 1.35% 16.83%

0.39%

44.34%

3.08E+10 1.22% 14.52%

0.35%

50.10%

3.23E+10 1.21% 21.58%

0.37%

0.74%

L1D
ミス率

L1D ミスdm 率

 チューニングにより、L2デマンドミス率が減少し、L2スループットが高い状態となった。
FUJITSU CONFIDENTIAL

12

Copyright 2018 FUJITSU LIMITEDGCRソルバのPost-FX100性能推定
 GCRカーネル
行列特性：非対称行列,17点ステンシル（5点差分✕4次元）
問題サイズ：NxxNyxNzxNv=48x48x16x128（16と12で割り切れるようにNxxNy=48x48と調整）
コンパイラオプション：(共通)-xpre̲diff -Kfast -Kopenmp -Kocl,optmsg=2 -V -X9 -Kmfunc=2 -Cpp -Kloop̲nofusion-CcdRR8 -Kswp̲strong -Kprefetch̲nostrong -Nlst=t,lst=i,lst=x Kprefetch̲sequential=soft(DAXPY) -Kprefetch̲cache̲level=all(SpMV) -Kprefetch̲cache̲level=1

FUJITSU CONFIDENTIAL

13

Copyright 2018 FUJITSU LIMITEDDAXPY1
 FX100 の実行時間と Post-FX100推定値の比較
レジスタ数

FX100 (PA実行時間)

128reg32reg

P-FX100 (推定値)

スレッド数

実行時間(sec)

比率

0.32

1.00

0.32

1.00

0.24

0.75

16
12

FX100
128reg
16t

FX100
32reg
16t

P-FX100
32reg 12t推定

 推定性能について メモリバンド幅ネックのループである。
Post-FX100推定：32reg 12t ASIS

浮動小数点ロードL2アクセス待ち

ビジー率

メモリビ メモリRビ メモリWビ L2ビジー L1ビジー
ジー率 ジー率 ジー率率率

Post-FX100 100%

71%

51%

87%

44%

浮動小数点ロード
メモリアクセス待ち

メモリ・キャッシュ
ビジー待ち

FUJITSU CONFIDENTIAL

14

ビジー待ち

Copyright 2018 FUJITSU LIMITEDDAXPY2
 FX100 の実行時間と Post-FX100推定値の比較
レジスタ数

FX100 (PA実行時間)
Post-FX100 (推定値)128reg
32reg

スレッド数

実行時間(sec)

比率

0.32

1.00

0.32

1.00

0.24

0.75

16
12

FX100
128reg
16t

FX100
32reg
16t

P-FX100
32reg 12t推定

 推定性能について
 VEC1 と同様に、メモリバンド幅ネックのループである。
浮動小数点ロードL2アクセス待ちPost-FX100推定：32reg 12t ASIS

ビジー率

メモリビ
ジー率

Post-FX100 100%

メモリRビ メモリWビ L2ビジー L1ビジー
ジー率
ジー率
率
率

72%

50%

87%

44%
浮動小数点ロード
メモリアクセス待ち

メモリ・キャッシュ
ビジー待ち

FUJITSU CONFIDENTIAL

15

ビジー待ち

Copyright 2018 FUJITSU LIMITEDSpMV FX100 の実行時間と Post-FX100 推定値の比較
コード

レジスタ数

asis

128reg

tune

32reg

FX100 (PA実行時間)
Post-FX100 (推定値)(2次元分割)

スレッド数

16
12

実行時間(sec)

比率

0.37

1.00

0.38

1.02

0.37

1.00

 推定性能について L2ビジー傾向のループとなっている。 128regと同等な性能となった。

FX100
128reg
16t

FX100
32reg
16t

P-FX100
32reg 12t

推定

Post-FX100 推定 ： 32reg-12t チューニング版

浮動小数点ロードL2アクセス待ち

ビジー率

メモリビジー メモリRビジー メモリWビ率
率
ジー率

Post-FX100 51%
FUJITSU CONFIDENTIAL

41%

18%

浮動小数点 整数演算パ アドレス計算
L2ビジー率 L1ビジー率 演算パイプラ イプラインビ パイプライン
インビジー率 ジー率
ビジー率

84%

70%

16

20%

11%

30%

メモリ・キャッシュ
ビジー待ち

Copyright 2018 FUJITSU LIMITED4.4 N 体カーネルの性能推定
理化学研究所計算科学研究センター 似鳥 啓吾
富士通株式会社次世代ＴＣ開発本部 青木 正樹
4.4.1 概要
前回の WG（ポストペタアプリ性能 WG、P-PAP）では倍精度演算でニュートン重力の時間導関数を含めた N 体計算を、時間積分を含めて性能評価を行った。対象としたアーキテク
チャは富士通 FX10/FX100、Intel KNC 等である。このコードでは


組み込み関数サポートを含むコンパイラ性能

演算数の少ないときのコア間同期性能を評価するようなベンチマークとなっている。今回はもう少し単純な、単精度で導関数を含まない重力のみのカーネルの評価を行った。このような条件にした理由は
SVE で単精度 SIMD がフルサポートされるようになったこと



レジスタ数が 32 本になった影響を見極めておきたかったこと



実機での評価はまだ先になりそうなのでできるだけ単純な条件にしておきたかったことが挙げられる。よって本 WG での評価では単一コアで命令列がいかに効率よく実行パイプラ
インを埋めることができるかに主眼を向けたものとなる。N 体計算のカーネルは主記憶帯域への要求は低いものの命令間の依存関係という点では DGEMM や SGEMM のような行列行列乗
算よりも性能を出すのが難しいケースもある。
4.4.2 Intel KNL での計測
第 2 回会合（2017/4/19）では、Intel KNL 上で AVX-512 命令を用いたコードの性能について発表した。よく知られているとおり、ポスト「京」での SIMD 幅も 512-bit であるため、コード移植の準備や比較対象のアーキテクチャとして KNL は適している。
この作業の中で最も驚きだったのが、Intel コンパイラが外側ループに対して#pragma
simd を付加しただけで理想的に SIMD 化されたコードを出してしまったことである。O(N2)の N 体計算カーネルには重力を受け取る粒子についての外側の i-ループと重力の源となる粒子についての内側の j-ループがある。通常コンパイラが自動ベクトル化を行うのは最内ループに対してであるが、計算効率の観点からは外側ループの並列度が利用可能なときはこちらを利用するのが望ましい。ディレクティブをひとつ付けただけでメモリオペランドの放送を含めてコンパイラが外側ループの SIMD 化をしてしまったこと、組み込み関数版も作ってはみたもののアセンブリ出力も性能もコンパイラベクトル化のものとほぼ同一であったことを書き留めておく。
KNL の命令セットでひとつ特徴的なのは、逆数や逆数平方根の近似命令に対して単一命令
で 28-bit 精度のものが用意されていることである。これらの命令は AVX-512ER というカテゴリに含まれ、Skylake Xeon から利用可能なのは 14-bit 精度の近似命令のみである。倍精度に対しては 1 回の Newton 法反復、単精度の場合は反復なしでフルの精度が得られる。追加の回路と引き換えに単に命令数が減った分だけ速くなるだけでなく、数値定数格納のためのレジスタや out-of-order リソースの消費削減も期待できる。
計測された性能であるが、KNL の 1 コアに対して 1/2/4 スレッドを用いた場合ループ回転のサイクル数にして 14.55/13.00/11.34 であった。ループ中の命令数が 18、内 AVX 命令が14 であったので、理想値は 2 命令発行のアーキテクチャの場合 9 サイクル、3 命令発行でレジスタ間の mov 命令を排除できたとすると 6.5 サイクルであり、良好なスケジューリン
グが得られているといえる。
少し別の話題にはなるが、この WG では SIMD レジスタを用いたテーブル参照の手法も紹介した。ARM SVE にも同様の tbl という命令が存在する。周期境界条件の分子動力学や宇宙論 N 体計算では相補誤差関数やガウス関数が相互作用カーネルに含まれるため、このような手法が威力を発揮する。
4.4.3 SVE 命令セットのコンパイラ出力
第６回会合（2018/4/16）では、富士通より SVE 命令用にコンパイルされた結果の提供を受け、その中身を簡単に紹介した。実機やシミュレータでの性能評価には至っていない。
前回の Intel コンパイラでは外側ループに対して SIMD 化が行われたが、富士通コンパイラにはそのような機能は未対応とのことで、内側の j-ループに対して SIMD 化とソフトウェア
パイプライニングが適用された。粒子データは 3 つの座標と 1 つの質量の 4 単精度語をパックした AoS (Array of Structures)になっていたため、粒子データのレジスタへの読み込
みが Gather load になってしまっている。外側ループの SIMD 化ではそれぞれの語をレジス
タないしメモリオペランドとして水平に放送していたので問題にならなかった箇所ではあるが、ここは SIMD 化の方針によっては SoA (Structure of Arrays)にユーザー側で変更すべき箇所であろう。逆数平方根の部分は、数値定数を 3 つ保持して有効桁数を一度に 3 倍にする方法が用いられていた。これは HPC-ACE のときの単精度逆数平方根の扱いと同様である。ARM SVE には数値定数のためにレジスタを消費せず有効桁数を 2 倍にする補助命令も用意されているがこちらは使われていなかった。
4.4.4 富士通での性能評価
第 8 回会合（2018/10/19）では富士通からの性能評価の報告があった。ただしこれは
A64FX の実機やシミュレータでの評価ではなく、FX100 でレジスタ数を 32 本に制限してコンパイル・実行したもの、もしくはそのプロファイリング情報から Post-FX100 での性能を推計したものとなる。単純にレジスタ数を削減した結果は 128 レジスタ時にはほぼ見られなかった浮動小数点演算待ちが大幅に増加し、性能は 49%まで低下してしまっている。これはレジスタ不足によりソフトウェアパイプライニングが適用されなくなったため（SVE は
swp されてる）で、その後ソフトウェアパイプライニングを適用すべく手動でループを 2 分割した版が示された。この結果ループ分割前から 22%の性能向上が得られ 128 レジスタ時の60%の性能となっている。
なお、単純に FX100 でレジスタ数を削減しただけでは A64FX での正確な性能予測は不可能であり、命令セットの変更によるアドレス計算の削減や強化された out-of-order リソー

スといった性能改善要素を加味する必要がある（現時点では直接数字を出すことはできな
いが、実際にシミュレータを用いた計測で推計値よりも性能が改善することが確認された）。
Out-of-order や緩やかなソフトウェアパイプライニングが性能改善をもたらした場合、
ループ分割の是非も逆転する可能性がある。この文書を執筆開始した 2019 年 2 月現在ではフラッグシップ２０２０プロジェクトにおいて DGEMM や STREAM といった重要なカーネルの実機評価の数値も徐々に公開され始めているので、遠くない将来に一般ユーザーにとっても実機のデータが入手可能になることを期待したい。

N体計算ベンチマーク（AVX-512で遊んでみた話）
2017/4/19
SS研メニーコア時代のアプリ性能検討WG第2回会合
理化学研究所計算科学研究機構
似鳥啓吾

はじめに
• 実機で試せる新アーキテクチャということでAVX-512 on KNL
を試してみた• 以下の２つをやってみた
1. 単純な重力多体計算カーネル• コードはhttps://github.com/nitadori/gravity-512/に

2. SIMDレジスタを用いたテーブル参照

AVX-512について• インテルのSIMD命令セットとしては大幅なメジャーアップ
デート• KNightsCornerでの拡張命令はx86 CPUのブランチ的な存在だったが、AVX-512ではSSE/AVXとの互換性も確保しつつ
Skylake Xeonからも採用• 現状ではKNightsLandingsで実機テストが可能• Skylake Xeon E5 v5はそろそろ発売？Amazon cloudで最初に触れるらしい• レジスタ幅が512-bitに、倍精度なら8語、単精度なら16語• レジスタの本数が32本へと倍増、マスクレジスタ新設• オペランド修飾が充実、特にメモリオペランドから放送ができる重力多体問題計算カーネル#pragma omp parallel for#pragma simdfor(int i=0; i<N; i++){float ax=0.0f;
float ay=0.0f;
float az=0.0f;
float po=0.0f;
float xi = posm[i].x;
float yi = posm[i].y;
float zi = posm[i].z;
for(int j=0; j<N; j++){float dx = xi ‐ posm[j].x;
float dy = yi ‐ posm[j].y;float dz = zi ‐ posm[j].z;
float r2 =eps2 + dx*dx + dy*dy + dz*dz;

float ri = 1.0f / sqrtf(r2);float mri= posm[j].m * ri;

float ri2

= ri * ri;

float mri3 = mri * ri2;
ax ‐= mri3 * dx;
ay ‐= mri3 * dy;
az ‐= mri3 * dz;}
accp[i].ax = ax;
accp[i].ay = ay;
accp[i].az = az;}}• 外側ループの並列度に対して#pragma simdを付けただけ• 内側ループの自動ベクトル化よりは高い性能が期待できる• 内側ループには逆数平方根

アセンブリ出力
1..B1.42:# Preds ..B1.42 ..B1.41

2# Execution count [2.50e+01]

3

vsubps(%rcx,%r10){1to16}, %zmm10, %zmm18

4

vmovaps%zmm4, %zmm14

5

vsubps4(%rcx,%r10){1to16}, %zmm9, %zmm19

6

addq

$1, %r11

7

vfmadd231ps %zmm18, %zmm18, %zmm14

8

vsubps9

vfmadd231ps %zmm19, %zmm19, %zmm14

10

vfmadd231ps %zmm21, %zmm21, %zmm14

11

vrsqrt28ps %zmm14, %zmm15{%k1}{z}12vmulps12(%rcx,%r10){1to16}, %zmm15, %zmm16

13

vmulps%zmm15, %zmm15, %zmm17

14

addq

$16, %r10

15vmulps%zmm17, %zmm16, %zmm20

16

vfnmadd231ps %zmm18, %zmm20, %zmm13

17

vfnmadd231ps %zmm19, %zmm20, %zmm12

18

vfnmadd231ps %zmm21, %zmm20, %zmm11

19

cmpq

%r15, %r11

20

jb..B1.428(%rcx,%r10){1to16}, %zmm8, %zmm21# Prob 82%• インテルコンパイラがほぼ完
璧なコードを吐いてしまった• 人生であまりない経験• 放送ロードがメモリオペラン
ドに入る• 単精度の逆数平方根は1命令• 18命令（内AVX命令は14）• 理想的にレイテンシが隠蔽されると9 cycle/loop程度か

組込み関数版• 書いては見たものの、#pragma simdによる自動ベ
クトル化版と全く変わらない
出力• この一件をもって「SIMD化
なんかコンパイラが全部やってくれる」と一般化された期待をするべきではないと思いますが、隔世の感はありますKNLでのシングルコア性能• 1.3 GHz 1コアでSMTを1/2/4にして計測• 1スレッドでもout-of-orderがうまく機能している• KNCでは最低2スレッドが必要でした• 4スレッドにすることで27.7%の性能向上

core/thread

nsec/loop

cycle/loop

1C1T

11.19

14.55

1C2T

10.00

13.00

1C4T

8.76

11.34

テーブル参照の話• 実数1入力1出力の数学関数高速化においては常套手段• 特に入力区間や必要精度がわかっているときに強力• 例：分子動力学での相補誤差関数や指数関数• 必ずしも最近のマシンで性能が出るわけではない• キャッシュに収めてなおかつヒットさせないとパイプラインがストー
ルする• SIMD化が実質不可能、gatherは遅い• Wide SIMDのレジスタに対する水平シャッフルでテーブル参照
を行ってみた• AVX-512だと単精度なら16エントリのテーブル• 2つのレジスタをソースに32エントリとできる命令まである• KNLで簡単なテストを作ってみた• float exp2f(float x), y=2x for x∈[1.0, 2.0)
• 3次の多項式

SIMDでテーブル参照• こんなことができる命令がある
b[i] = a[idx[i]];• aとbはfloat[16]のSIMDレジスタ• idxはint[16]のレジスタ(下４-bitのみを使う 0<=idx[i]<16)• AVX-512の該当命令は：vpermps、組込み関数だと__m512 _mm512_permutexvar_ps (__m512i idx, __m512 a)

もっと強力な命令もある• vpermi2ps, vpermt2p• __m512 _mm512_permutex2var_ps (__m512 a, __m512i idx, __m512 b)• データにレジスタ２本、インデクスに１本• aとbを繋いでidxで要素を参照
c[i] = idx&16 ? b[idx[i&15]] : a[idx[i&15]];• 単精度32エントリのテーブルとして利用可能• FMA3同様、３つの入力オペランドのうち１つが破壊される

Cではこんなコード• float exp2f(float x), y=2x for x∈[1.0, 2.0)
• 区間を32分割して3次の多項式で近似• 単精度の丸め誤差程度までの精度が得られた• 結果的に16エントリでも十分だった
float exp2app_32(const float x){
static const float tab[32][4] = {#include"tab32x4.h"};
union{
floatf;unsigned u;} m32;
m32.f = x;m32.u >>= 18;
const int idx = m32.u & 31;
m32.u <<= 18;
const float h = x - m32.f;
const float *a = tab[idx];
const float y = a[0] + h * (a[1] + h * (a[2] + h * (a[3])));
return y;}テーブルの生成• SollyaというツールのRemezアルゴリズム(minimax近似)を用いた• http://sollya.gforge.inria.fr/
prec = 60;
f = 2^x;
for n from 0 to 31 do {p = remez(f(x+1+n/32), 3, [0;1/32]);// p;
print("{",
coeff(p,0), ",",
coeff(p,1), ",",
coeff(p,2), ",",
coeff(p,3), "}," );print("//", dirtyinfnorm(p - f(x+1+n/32), [0;1/32]));};

実装• さすがにこのぐらい複雑なものになると組込み関数が必要• 原理的には特定の大きさのテーブル参照に対してはコンパイラからの生成も可能
for(int i=0; i<n; i+=16){__m512 xv = _mm512_loadu_ps(&x[i]);__m512i sr = _mm512_srli_epi32((__m512i)xv, 18);__m512i sl = _mm512_slli_epi32(sr, 18);__m512 dx = _mm512_sub_ps(xv, (__m512)sl);//__m512i idx = _mm512_and_epi32(sr, _mm512_set1_epi32(31));__m512i idx = sr;__m512 a0 = _mm512_permutex2var_ps(z00, idx, z01);__m512 a1 = _mm512_permutex2var_ps(z10, idx, z11);__m512 a2 = _mm512_permutex2var_ps(z20, idx, z21);__m512 a3 = _mm512_permutex2var_ps(z30, idx, z31);__m512 yv =_mm512_fmadd_ps(dx,_mm512_fmadd_ps(dx,_mm512_fmadd_ps(dx, a3, a2), a1), a0);_mm512_storeu_ps(&y[i], yv);}

vmovups(%r14,%rax,4), %zmm8
vmovaps%zmm7, %zmm13
vmovaps%zmm5, %zmm12
vmovaps%zmm3, %zmm11
vpsrld$18, %zmm8, %zmm10
vpslld$18, %zmm10, %zmm9
vpermt2ps %zmm6, %zmm10, %zmm13
vsubps%zmm9, %zmm8, %zmm14
vpermt2ps %zmm4, %zmm10, %zmm12
vpermt2ps %zmm2, %zmm10, %zmm11
vpermi2ps %zmm0, %zmm1, %zmm10
vfmadd231ps %zmm14, %zmm10, %zmm11
vfmadd231ps %zmm14, %zmm11, %zmm12
vfmadd213ps %zmm13, %zmm12, %zmm14
vmovups%zmm14, (%rbx,%rax,4)
addq$16, %rax
cmpq%r12, %raxjl..B1.30

性能• On 1 core of Intel(R) Xeon Phi(TM) CPU 7210 @ 1.30GHz• 32768 elements, 256 KiB read and 256 KiB write• 16-entry, 1-thread• 0.487948 nsec/element, 7.807161 nsec/loop• 16-entry, 4-thread• 0.424348 nsec/element, 6.789574 nsec/loop• 32-entry, 1-thread• 0.661530 nsec/element, 10.584481 nsec/loop• 32-entry, 4-thread• 0.634580 nsec/element, 10.153279 nsec/loop

まとめ• SIMD幅がどんどん広がると使い勝手が悪化していくというのが一般論だけど• 広大なレジスタ空間（zmmが32本あると2 KiB、4スレッドで8KiB）を活用するという楽しみ方もあります• とはいっても、Xeon PhiではなくてSkylake XeonのAVX-512の方が本命だと思います重力多体問題（FujitsuコンパイラSVE出力）
2018/4/16
SS研メニーコアWG
理研CCS 似鳥啓吾

概要
• 富士通様から重力多体問題計算カーネルのSVE向けコンパイル
結果を提供していただきました• 動かして性能計測はまだですが、コンパイル結果を見ていこうと思いますARM SVEについて• Scalable Vector Extension• ARM社のSIMD拡張命令• 特徴：• Vector Length Agnostic (VLA)：ベクトル長不可知論• 128-bitの倍数のSIMD幅なら最大2048-bitまでSIMD幅は「実装依存」• 同じバイナリが複数の異なるチップで動く• SIMD長の実装に依りループの反復数が変化する• （適切にコンパイルされていれば）結果はSIMD長に依らず一致• Intelと比較して• 同じAVX-512でもKNLとSkylakeでサポートしている命令が違う• ただしKNxは市場から消え去りそう• Fujitsu SPARC HPC-ACEと比較して• 倍精度にあったあの命令が単精度には無い、みたいなことが起こらない• すべての命令を検証するのは大変そう・・・入力C++コード• 普通の2重ループ• ただしデータ構造はAoS• 外側ループに対するSIMD化と放送ロードの使用を想定しました• 内側ループに対してSIMD化す
るとstructure loadになりますコンパイラのメッセージ• Optimization messages•

jwd6004s-i "test1.cpp", line 34: リダクション演算を含むループ制御変数'j'のループをSIMD化しました。

•

jwd8204o-i "test1.cpp", line 34: ループにソフトウェアパイプライニングを適用しました。

•

jwd8205o-i "test1.cpp", line 34: ループの回転数が96回以上の時、ソフトウェアパイプライニングを適用したループが実行時に選択されます。

•

jwd8209o-i "test1.cpp", line 39: 多項式の演算順序を変更しました。

•

jwd8206o-i "test1.cpp", line 40: 除算を逆数の乗算に変更しました。

• 内側ループに対してSIMD化とソフトウェアパイプライ
ニングが行われましたコンパイル結果• SWPがかかったものはちょっととても読めない
ので、端数用の単純ルー
プ• SWPの方もspillはなかった• 残念ながらGather load
命令になっていまっている• ld1wという命令• まぁこれは元のデータ構造の側の問題でもある• movprfxについては後述

SVEでのFMA (fused multiply-add)• 命令フィールドの都合からか3 operandだけ取る• 加減算もしくは乗算オペランドのどちらかが破壊される• Intel AVX2のFMA3と同様• 加減算オペランドのレジスタを破壊する命令• FMLA, FMLS, FNMLA, FNMLS• 乗算オペランドのレジスタを破壊する命令• FMAD, FMSB, FNMAD, FNMSB• どちらのオペランドも破壊されて欲しくないとき• movprfx命令を直前に書いて破壊されるレジスタを避難させる• ただし実行コストを伴わない• 後続命令を修飾して4オペランドFMAのように振る舞わせる• こういった知識は通常のアプリプログラミングでは不要。コン
パイラの出力を読むためにはある程度慣れておく必要がある。

逆数平方根の部分• 有効桁数を3倍にする公式が使われている• 初期値推定の精度は8-bit程度であると予想される• この点はHPC-ACEと同じ• SVEにはRSQRTSという収束補助の命令がある• 1/2倍するのに余計な命令を消費しない• 3.0のような数値定数の保存にレ
ジスタを消費しない• 倍精度のときはこれを3回呼ぶのだと思う逆数平方根命令略史• AMD 3DNow!（20周年）• float2 SIMDに対して15-bit精度• 収束補助命令が2つ• 21個の新命令に対して2が逆数、3つが逆数平方根まわりだったことになる• Intel SSE• float4 SIMDに対して12-bit精度• このあとず〜〜っとAVX-512になるまで倍精度からは近似命令を呼べなかった• IntelとAMDでビットレベルの互換性はなかった• HPC-ACE• double2もしくはfloat2に対して8-bit精度• HPC-ACE2• double4もしくはfloat4。float8は？

• AVX-512• 14-bit精度になった、倍精度からも単精度からも利用可能• KNLではいきなり28-bit精度が得られる拡張が利用可能

要望• 外側ループに対してSIMD化したい• AoSのまま放送ロードで対応できる• VLAの思想的にも、総和を含むSIMD化はベクトル長ごとに結果が異なってしまう• できればディレクティブで• Intelコンパイラにはでたことなので• OpenMP 4.0かOCLかはどちらでもいいけど• 今のコンパイラで既に実現可能ならコードを見てみたい• まだできないなら実機を我々が触れる時期までに• これができるようになるとFDPS on Post-Kが素晴らしく使いやすいものになりますSS研メニーコア時代のアプリ性能検討WG資料

N体 カーネルの評価報告

2018年10月19日
富士通株式会社
Copyright 2018 FUJITSU LIMITED目次
 測定環境 評価環境 FX100 の測定条件と評価コード

 FX100 の性能分析
 最適化の状況
 128reg と 32reg コンパイラ の比較
 チューニング Post-FX100 推定 Post-FX100 のアーキテクチャ
 性能予測ツールの概要
 評価区間の推定 (チューニング)

 まとめ

1

Copyright 2018 FUJITSU LIMITED

測定環境
 評価環境 FX100 の測定条件と評価コード2

Copyright 2018 FUJITSU LIMITED評価環境 FX100 (ハードウェア、ソフトウェアの情報)
FX100 諸元

C
P
U

名称
動作周波数
コア数

SPARC64 Xifx
1.975GHz
32コア＋アシスタントコア

キャッシュメモリ容量

L1I$ : 64KB/コア(4way)
L1D$: 64KB/コア(4way)
L2$ : 24MB/CPU

理論ピーク性能
主記憶-CPU間帯域
CPU数、コア数ノー
ド理論ピーク性能

1011GFlops/CPU(倍精度)2022GFlops/CPU(単精度)(240+240)GB/s/CPU
1CPU/ノード(16コアｘ2CMG＋アシスタントコア)/ノード
1011GFlops/ノード(倍精度)2022GFlops/ノード(単精度)

主記憶容量

32GB/ノード

主記憶-CPU間帯域
ノード間通信帯域
OS の名称とバージョン(240+240)GB/s
100GB/s
FX100向けOS 2.1.1
3

Copyright 2018 FUJITSU LIMITEDFX100 の測定条件と評価コード
 FX100 の測定条件
項目

条件

コンパイラ

Post-FX100向けチューニング用コンパイラ

翻訳時オプション

ASIS ： -Kfast,ocl,openmp 、ループ分割時 ： -Kloop_nofusion

スレッド数

16スレッド

実行条件

1ノード – 1プロセス

PA 採取条件

評価区間をイタレーション1,000で回転させてPA情報を採取

 評価コード 似鳥様からご紹介頂いた
https://github.com/nitadori/gravity-512
よりダウンロードしたソースを利用 評価区間は gravity_simple
 規模 N は 1,024 (オリジナル ： N = 16 * 1,024) 単精度 (C++)

4

Copyright 2018 FUJITSU LIMITEDFX100 の性能分析
 最適化の状況
 128reg と 32reg コンパイラ の比較
 チューニング

5

Copyright 2018 FUJITSU LIMITED

最適化の状況 (評価区間のソースリスト) ソースコードの概要 (gravity_simple)
:
26
27
28
29
30
31
32
33
34 p
35 p
36 p
37 p
38
39
40 p
41 p
42 p
43

void gravity_simple (const int N,
const float eps2,
const PosM * __restrict posm,
AccP * __restrict accp){#pragma omp parallel for// #pragma simd
(現状 omp simd は未対応)
for (int i=0; i<N; i++) {float ax=0.0f;
float ay=0.0f;
float az=0.0f;
float po=0.0f;
float xi = posm[i].x;
float yi = posm[i].y;
float zi = posm[i].z;<<< Loop-information Start >>><<< [OPTIMIZATION]
ループ回転数128回以上で適用<<< SIMD(VL: 4)<<< SOFTWARE PIPELINING<<< Loop-information End >>>44 p 4v for (int j=0; j<N; j++) {45 p 4v
float dx = xi - posm[j].x;
ループ回転数N は 1024
46 p 4float dy = yi - posm[j].y;
47 p 4float dz = zi - posm[j].z;
48
49 p 4float r2 = eps2 + dx*dx + dy*dy + dz*dz;
50 p 4v
float ri = 1.0f / sqrtf(r2);5152 p 4v
float mri = posm[j].m * ri;
53 p 4v
float ri2 = ri * ri;
54 p 4v
float mri3 = mri * ri2;
55
56 p 4v
ax -= mri3 * dx;
57 p 4v
ay -= mri3 * dy;
58 p 4v
az -= mri3 * dz;
59 p 4v }60
61 p
accp[i].ax = ax;
62 p
accp[i].ay = ay;
63 p
accp[i].az = az;
64 p
}
65
}
6

Copyright 2018 FUJITSU LIMITED128reg と 32reg コンパイラ の比較
 PAの実行時間で比較
 128reg と 32reg コンパイラの性能向上比は 0.49倍
であり、レジスタ数減少による性能影響は大きい。
gravity_simple

レジスタ数

スレッド数

128

実行時間
(sec)

性能
向上比

0.153

1.00

0.313

0.49

16

FX100

32

FX100
128reg
16t

FX100
32reg
16t

4 命令コミット

2/3 命令コミット

1 命令コミット

 最適化状況 32reg では、レジスタ不足でソフトウェアパイプライニングが適用されないため、浮動小数点演算待ち が増加し、性能に影響が出ている。 ⇒ 32reg ではループ分割が必要。
レジスタ数

SIMD

UNROLL

swpl 適用

Reg Fp128

44128回以上85 / 128

32

4

2

浮動小数点
レジスタ不足

－

7

浮動小数点
演算待ち

Copyright 2018 FUJITSU LIMITEDチューニング ： ループ分割の検討
 ループの依存関係図から分割点を検討
最内ループ：LOOP_L217

２

３

１

for 文の最内ループ(L217) では、１ 、 ２ 、 ３ の位置で、 (図中の線で表す)レジスタへの依存が少ない。

Post-FX100 のアセンブラ
をツールを用いて図化
33 void gravity_simple (const float eps2)
34 {35 #pragma omp parallel for
36 // #pragma simd
外側ループ：L214
37
for (int i=0; i<N; i++) {38
float ax=0.0f;
39
float ay=0.0f;
40
float az=0.0f;
41
float po=0.0f;
42
43
float xi = posm_x[i];44float yi = posm_y[i];45float zi = posm_z[i];46
最内ループ：L217
47
for (int j=0; j<N; j++) {48
float dx = xi - posm_x[j];49
float dy = yi - posm_y[j];50float dz = zi - posm_z[j];
51

52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
}
68 }

float r2 = eps2 + dx*dx + dy*dy + dz*dz;⇒ 候補 １
float ri = 1.0f / sqrtf(r2);⇒ 候補 ２
float mri = posm_m[j] * ri;
float ri2 = ri * ri;
float mri3 = mri * ri2;
⇒ 候補 ３ax -= mri3 * dx;
ay -= mri3 * dy;
az -= mri3 * dz;}
accp_ax[i] = ax;
accp_ay[i] = ay;
accp_az[i] = az;8

以下の観点から、候補2 で評価を実施する・使用レジスタ数が少ない
・アンロールが必要となるリダクション部分と使用レジスタ数が多い 1/sqrtf部分 を分離
・命令数バランスを考慮
Copyright 2018 FUJITSU LIMITEDチューニング ： ソース (最内ループを2分割)
 チューニングのソースコードの概要 (gravity_simple)
:
38
39
40 p
41 p
42 p
43 p
44
45 p
46 p
47 p
48
49

59

// loop-div-2<<< Loop-information Start >>><<< [OPTIMIZATION]
ループ回転数48回以上で適用<<< SIMD(VL: 4)<<< SOFTWARE PIPELINING<<< Loop-information End >>>60 p 3v for (int j=0 ; j<N ; j++) {61 p 3v
float ri = wk_ri[j];ループ回転数62
N は 1024
63 p 3float mri = posm[j].m * ri;
64 p 3v
float ri2 = ri * ri;
65 p 3v
float mri3 = mri * ri2;
区間２
66
67
// work (x,y,z)
68 p 3v
ax -= mri3 * wk_dx[j];69 p 3v
ay -= mri3 * wk_dy[j];70 p 3v
az -= mri3 * wk_dz[j];71 p 3v }72
73 p
accp[i].ax = ax;
74 p
accp[i].ay = ay;
75 p
accp[i].az = az;
76 p
}
77
}

追加した work 配列#pragma omp parallel for
private(wk_dx,wk_dy,wk_dz,wk_ri)// #pragma simd
(現状 omp simd は未対応)
for (int i=0; i<N; i++) {float ax=0.0f;
ループ回転数 N は 1024
float ay=0.0f;
float az=0.0f;
float xi = posm[i].x;
float yi = posm[i].y;
float zi = posm[i].z;// loop-div-1<<< Loop-information Start >>><<< [OPTIMIZATION]
ループ回転数32回以上で適用<<< SIMD(VL: 4)<<< SOFTWARE PIPELINING<<< Loop-information End >>>50 p 2v for (int j=0 ; j<N ; j++) {51// work (x,y,z)
ループ回転数52 p 2v
wk_dx[j] = xi - posm[j].x;
N は 1024
53 p 2v
wk_dy[j] = yi - posm[j].y;
54 p 2v
wk_dz[j] = zi - posm[j].z;
区間１55
56 p 2
float r2 = eps2 + wk_dx[j]*wk_dx[j]+ wk_dy[j]*wk_dy[j] + wk_dz[j]*wk_dz[j];57 p 2v
wk_ri[j] = 1.0f / sqrtf (r2);58 p 2v }赤字は ASIS からの修正点

9

Copyright 2018 FUJITSU LIMITEDコードチューニング
 PAの実行時間で比較
 ループ分割のチューニングにより、32reg ASIS とループ分割を比較すると 1.22倍 の性能向上が見られる。
gravity_simple

スレッド
数

レジスタ
数

ソース

128
FX100

16

32

ASIS
ﾙｰﾌﾟ分割

実行時間
(sec)

性能
向上比

0.153
0.313
0.255

1.00
1.22

ASIS
FX100
128reg16t

 最適化状況

1 命令コミット

ASIS
FX100
32reg
16t

tune
FX100
32reg
16t

4 命令コミット

 ループ分割によりソフトウェアパイプライニングが適用されたことで、浮動小数点演算待ちが減少して性能改善した。作業
領域の受け渡しや、アドレス計算などが増加したため、命令数は増加している。
2/3 命令コミット

ソース

ASIS
ﾙｰﾌﾟ
分割

レジスタ数

128
32

ループ分割

浮動小数点
演算待ち

SIMD

UNROLL

swpl 適用

Reg Fp4

4128回以上85 / 128

42

レジスタ不足

－

区間１

4

2

32回以上

30 / 32

区間２

4

3

48回以上

29 / 32

－

10

浮動小数点
ロード L1D
アクセス待ちCopyright 2018 FUJITSU LIMITEDPost-FX100 推定 Post-FX100 のアーキテクチャ
 性能予測ツールの概要
 評価区間の推定 (チューニング)

11

Copyright 2018 FUJITSU LIMITEDPost-FX100 のアーキテクチャ

SIMD

機能

Post-FX100

PRIMEHPC
FX100

京

ISA

ARMv8-A+ SVE

SPARC V9
+ HPC-ACE2

SPARC V9
+ HPC-ACE

bit幅

512

256

128

単精度要素数16

82

倍精度要素数842（命令に制限有）汎用整数レジスタ（本）

31+SP

32+32

32+32

（v8-Aと同じ）（g,l,i,o+xg）（g,l,i,o+xg）

ベクトルレジスタ（本）

32

128

128

プレディケートレジスタ（本）16（ベクトルレジスタを共用）（ベクトルレジスタ
を共用）

12

12

Copyright 2018 FUJITSU LIMITED

性能予測ツールの概要
アプリコード（カーネル部）
FX100：32regコンパイラを使用し詳細プロファイラ実行でPA情報を取得

PA情報：詳細プロファイラの出力結果(csvファイル）
性能予測ツール
性能予測値
 性能予測ツール概略
 FX100のPA情報を入力とする
 ハード、命令セット／コンパイラによる性能向上効果を算出する
 性能予測値を出力する

13

Copyright 2018 FUJITSU LIMITED

評価区間の推定 (コードチューニング)
 Post-FX100 推定値の比較 (ASIS と ループ分割)
gravity_simple

ソース

Post-FX100(単精度推定値)

ASIS
ﾙｰﾌﾟ分割

レジスタ
数

スレッド
数

32

12

実行時間
(sec)

性能
向上比

0.165

1.00

0.151

1.09

 推定性能について
 ループ分割を実施したことにより、FX100の結果同様ASIS性能に対して性能が改善した。

ASIS
Post-FX100
32reg
12t (推定) FX100の性能向上比と比較して、Post-FX100の向上比が小さい原因は、ループ分割による作業配列アクセス命令（アドレス計算）が増加し、命令コミットが増加したためである。
 ISA対応のシミュレーションでチューニングコードを性能測定した結果、チュー
ニング性能（推定）に対して約20%の性能向上を確認した。

tune
Post-FX100
32reg
12t (推定)4 命令コミット
2/3 命令コミット

削減

1 命令コミット

浮動小数点
演算待ち

削
減

Post-FX100実機では、更に高速化される可能性あり

14

Copyright 2018 FUJITSU LIMITED

まとめ

15

Copyright 2018 FUJITSU LIMITEDまとめ ASIS の問題点
 演算ネックのループのため、アーキレジスタ数減少の影響は出ている。
• 32reg コンパイラでは、レジスタ不足により、ソフトウェアパイプライニングが適用されず、浮動小数点演算待ちが増加し、性能に影響が出る。

 チューニングによる対応
 32reg コンパイラでは、ループ分割が必要となる。
• Post-FX100 では、コンパイラの分割機能により、自動で分割される。（予定） ループ分割によりソフトウェアパイプライニングが適用されて性能は改善した。
• ISA対応のシミュレーションでチューニングコードを性能測定した結果、チューニング
性能（推定）に対して約20%の性能向上を確認できた。

16

Copyright 2018 FUJITSU LIMITEDCopyright 2018 FUJITSU LIMITED