第2回 配信講義　計算科学技術特論B （2024）

計算科学技術 特論Ｂ 第2回 アプリケーションの性能最適化１ (高並列性能最適化) 2024年4月 ジャパンメディカルデバイス 開発本部�本部長 南�一生 [email protected] 2024年4月 計算科学技術特論B 講義全体の概要 • スーパーコンピュータとアプリケーションの性能 • アプリケーションの性能最適化1(高並列性能最適化) • アプリケーションの性能最適化の実例1とCPU単体性能と は？ • アプリケーションの性能最適化２(CPU単体性能最適化) • アプリケーションの性能最適化の実例2 2024年4月 計算科学技術特論B 2

今回の講義内容 • • • • 反復法(離散化・アルゴリズム) アプリケーションの性能最適化(並列化) 並列性能のボトルネック(並列化) 簡単な実例(並列化) 2024年4月 計算科学技術特論B 3 スーパーコンピュータを使うために (1)定式化(物理・数学の世界) (2)離散化(+-×÷の世界) (3)アルゴリズム(数学の世界) (4)プログラミング(コンピュータの世界) (5)並列化 (6)CPU単体性能チューニング 2024年4月 2022年4月 計算科学技術特論B 4

✦前回の講義で「高並列化」と「CPU単体性能の 向上」が重要なことを示しました ✦その両方について「並列処理」を利用しない と性能向上しない ✦「並列処理」できるアルゴリズムを採用する ことが重要とお伝えしました ✦最初の話題は「並列処理」できるアルゴリズ ムの話題です 2024年4月 計算科学技術特論B 5 ２次元偏微分方程式-離散化 左図のような２次元の領域で微分方程式を解くと する． ５点差分で差分化すると自分自身の差分点と周 りの4つの差分点の式ができる． C1,3u1 + C1,4 u2 + C1,5u7 = α 1 ①の差分点について C2,2u1 + C2,3u2 + C2,4 u3 + C2,5u8 = α 2 ②の差分点について ・・・・・・・・・・ C16,1u10 + C16,2u15 + C16,3u16 + C16,4 u17 + C16,5u22 = α 16 ・・・・・・・・・・ 行列の形 (疎行列) 2024年4月 計算科学技術特論B 6 ⑯の差分点について

２次元偏微分方程式-離散化 Cij 2024年4月 計算科学技術特論B 離散化 (+-×÷の世界) 7 離散化-離散化と連立一次方程式 左図のような２次元の領域で微分方程式を解くと する． ５点差分で差分化すると自分自身の差分点と周 りの4つの差分点の式ができる． C1,3u1 + C1,4 u2 + C1,5u7 = α 1 ①の差分点について C2,2u1 + C2,3u2 + C2,4 u3 + C2,5u8 = α 2 ②の差分点について ・・・・・・・・・・ C16,1u10 + C16,2u15 + C16,3u16 + C16,4 u17 + C16,5u22 = α 16 ⑯の差分点について ・・・・・・・・・・ 全部で42点の差分点が存在するので42元の連立方程式となる． したがって元の微分方程式は，離散化により以下のような連立一 次方程式を解く事に帰着させる事ができる． Ax = b 2024年4月 計算科学技術特論B 8

反復法 2024年4月 計算科学技術特論B 9 反復法の一般的な議論 Ax = b A = I−B 左のような連立一次方程式を考える． Aを右のように書き直す．I は単位行列． ・・・② (I − B)x = b ②を①に代入し整理する． x = Bx + b ③を満たすベクトルの列： を考える． ・・・① x 0 ,x1 ,x 2 ,x 3 ・・・ ・・・③ x (m+1) = Bx (m ) + b ・・・④ ここで④の両辺の極限を取れば⑤のようになる lim m→∝ x (m+1) = Blim m→∝ x (m ) + b・・・⑤ 極限ベクトルをxとおけば⑤は⑥の様になる x = Bx + b ・・・⑥ ⑥は③と全く同一となる．つまり連立一次方程式：①を解く事は④を 満たすベクトルの列： x 0 ,x1 ,x 2 ,x 3 ・・・ を見つける事に帰着する． 2024年4月 計算科学技術特論B 10

各反復法 x (m+1) = Bx (m ) + b ・・・④ ④を解くための反復法について説明する Aは正則，Aの対角要素は全て0でない数であると仮定する． 解くべき方程式 Ax = k ・・・⑦ Aを⑧のように変形．DはAの対角成 分のみの行列．EはAの狭義下三角 行列の符号を変えた行列．FはAの狭 義上三角行列の符号を変えた行列． A = D−E−F ・・・⑧ ⑧を⑦に代入し⑨を得る Dx = (E + F)x + k ・・・⑨ x (m+1) = D −1 (E + F)x (m ) + D −1k ⑨を④の反復法にすると⑩を得る． n この反復法を点ヤコビ法という． (m+1) ) aii xi = − ∑ aij x (m + kj 具体的には⑪の計算法となる． j j=1 注目すべきは(m+1)世代の解を得 j≠i るために(m)世代の解しか使用しな い点である． n ⎛ aij ⎞ (m ) ki (m+1) xi ・・・⑩ ・・・⑪ = −∑ ⎜ ⎟ x j + aii j=1 ⎝ aii ⎠ j≠i 2024年4月 計算科学技術特論B 11 各反復法 (D − E)x = Fx + k ・・・⑫ (D − E)x (m+1) = Fx (m ) + k ・・・⑬ ⑨を⑫のように変形する． ⑫を④の反復法にすると⑬を得る． この反復法を点ガウス・ザイデル法という． 具体的には⑭の計算法となる． i−1 (m+1) ii i ax = − ∑ aij x j=1 2024年4月 計算科学技術特論B 12 (m+1) j − n ∑a x ij j=1+1 (m ) j + kj ・・・⑭

うという欠点があります。 点ガウス・ザイデル法の参照関係 一方で、陰解法は未来の状態に基づく方程 このページに「いい i−1 n 式を考えて、解を求める方法です。計算のイ (m+1) (m+1) (m ) a x = − a x − a x + kj ii i ij j ij j メージを示すと図6.43のようになります。 j=1 j=1+1 株式会社ソ (m+1) 同じ(m+1)世代の別 レイドル ∑ ∑ A x の点を参照する． 約11ヶ月前 (m+1)世代 iの下位側の点を参照 (m)世代 帯行列 iの上位側の点を参照 ループ回転方向 図6.43 陰解法の計算イメージ (例)右図のように, iを定義するのにi-1の値を参 照するコーディングは下図のような参照関係とな り,リカレンス(回帰参照)が発生している． do i=1,n x(i) = a(i)* x(i − 1) + b(i) これを式の形で表現すると以下のようにな リカレンスがあると, iについて は依存関係があり並列処理が ります。 MESSENGER できない． 2024年4月 計算科学技術特論B 株式会社ソフトウ 13 リカレンスがあると．．． コメント 5 性能向上 大昔の計算機 (1)そのまま実行すると６clock 演算器 演算1 この計算の時点では 未来 の値は未知数で (2)演算を3つのス テージに分割する 演算2 演算i 株式会社ソ 演算1 す。したがって、未知数が複数含まれている (3)並列に計算できるようスケジューリング (4)２個の演算資源を使い並列に実行 演算3 演算2 この式をこのまま解くことはできません。そ (5)４clockで実行可 れでは、どのようにすれば解けるのでしょう 演算器に逐次に入力し逐次実行 か。もう少し具体的に見ていきたいと思いま (6)演算を細かく分割し複数の演算資源を 並列に動作し性能アップ す。 演算[i+1] 例として、図6.44に示すように1次元の領 処理のリカレンスをなくすチューニングは 域を5分割し、左から順に1, 2, 3, 4, 5と番号 現代のコンピュータを使う上で必須事項! をつけて考えます。 2024年4月 計算科学技術特論B 14 演算[i]

反復法のリカレンスの解消 ヤコビ法やガウス・ザイデル法はメッシュの計算順序を変更しても収束する事が数学的 に証明されている． ガウス・ザイデル法オリジナルの計算順序は前述のようにリカレンスのために並列化で きない． そこで(b)のハイバープレーン法や(c)のRED-BLACK法を用いループ内のリカレンス による依存関係をなくす。 依存関係をなくす事で前回講義のようにハードウェアの並列計算機構を使い高速に 計算可能となる．またコア間の並列化そのものにも使用可能． (a)オリジナルスイープ 参照関係 計算順序(スイープ順序) 2024年4月 計算科学技術特論B (b)ハイパープレーンスイープ 斜めのハイパープレーン内は リカレンスがない． ハイバープレーン間のみにリ カレンスが発生． (c)RED-BLACKスイープ 黄と青の２つのループに再構 成する． それぞれのループ内はリカレ ンスがない． 15 反復法に関する参考書 「計算機による大型行列の反復解法」 Ｒ．Ｓ．バーガ／著 渋谷政昭／〔ほか〕訳 2024年4月 計算科学技術特論B 16

共役勾配法(CG法) ①の2次の関数を考える． この関数の最小値問題を解くためにf(x)の 微分を取り0とおく． ②より①の最小値問題を解く事がax=bの解 を得る事と同等であることが分かる． 1 2 ax − bx 2 ・・・① f '(x) = ax − b = 0 ・・・② f (x) = 1 上記と同様な事を多次元のベクトルで行う． f (x) = (x,Ax) − (b,x) ・・・③ 2 ①と同等な関数として③を考える． n ③を成分で標記すると④のようになる． 1 n n f (x) = ∑ ∑ aij xi x j − ∑ bi xi ・・④ (Aは実対称行列であるとする) 2 i=1 j=1 i=1 ∂ f (x) 1 n 1 n = ∑ aik xi + ∑ akj x j − b・・⑤ k ∂xk 2 i=1 2 j=1 ④をxで偏微分すると⑤になる． 行列Aの対称性を使い０とおくことにより⑥ が得られる． ∂ f (x) n ⑥は連立一次方程式⑦の成分表示である． ∂x = ∑ aki xi − bk = 0 k i=1 Ax = b したがって③の最小値問題を解く事が連立 一次方程式⑦を解くこととなる． 2024年4月 計算科学技術特論B 17 共役勾配法(CG法) 1 f (x) = (x,Ax) − (b,x) 2 ・・・③ ③を求めるためのアルゴリズムについては 例えば ＜シリーズ新しい応用の数学 17＞ 「共役勾配法」 戸川隼人著 2024年4月 計算科学技術特論B 18 ・・・⑥ ・・・⑦

共役勾配法(CG法) CG法アルゴリズムの基本形を左図に示 す． この他に色々なアルゴリズムの派生形 がある． ③式のf(x)が最小値が得られるように xiの列を求めて行く． ri は残差ベクトルであり互いに直交し一 次独立であることが証明されている． したがって残差ベクトルはN元の連立一 次方程式にはN個しか存在しない． CG法を使うとN元の連立一次方程式は 高々N回の反復で収束する． またAの固有値が縮重しているか密集 していれば収束が速くなる性質を持つ． 2024年4月 計算科学技術特論B 19 前処理付き共役勾配法(CG法) Aに近い正値対称行列:Mを考えコレスキー分 解する（①式）． ここで②のようなxの置き換えを行う． M = UT U x! = Ux ・・・① ・・・② これを使ってもとのAx=bを変形する（③④式） AU −1x! = b ・・・③ -T (Ax=bに②を代入すれば③となる．③にU を 左辺からかければ④になる) U −T AU −1x! = U −T b = b! ・・④ ここで⑤で置き換えたAの性質をみてみる． MはAに近いので⑥を満たす． したがって①より⑦を満たす． ! U −T AU −1 = A A≈M A ≈ UT U ・・⑤ ・・・⑥ ・・・⑦ ⑤に⑦を適用すると⑧をみたし置き換えたA U −T AU −1 ≈ U −T (U T U)U −1 = I の性質は単位行列に近いことが分かる． ・・・⑧ したがって⑨によって置き換えた行列Aの固 有値は1の周りに密集しているものと考えら れる 置き換えた行列とベクトルに対する連立一次 方程式：⑩の収束は早いものと期待できる． 2024年4月 計算科学技術特論B 20 ! = U −T AU −1 A ! x! = b! A ・・⑨ ・・・⑩

前処理付き共役勾配法(CG法) Aは例えば先例に示した５点差分で得られる疎行列である とする． ここでAのコレスキー分解を考える． Aを完全にコレスキー分解すると方程式を解くことになる． Aの完全なコレスキー分解を行った行列は密行列となる． Aの完全なコレスキー分解のための演算量は大きい． そこでAを完全でなく不完全なコレスキー分解を行うこととす る． T A≈U U 不完全コレスキー分解とは元の行列:Aが持つ要素のみ出現 するようにコレスキー分解を行う． このような不完全コレスキー分解を実施することを前処理と 呼ぶ． 前処理した行列について共役勾配法を使用して連立一次 方程式を解く方法を前処理付き共役勾配法という． Aの完全コレスキー分解 Aの不完全コレスキー分解 Ã x̃ = b̃ Ã = U −T AU −1 2024年4月 行列:A 計算科学技術特論B 21 前処理付き共役勾配法(CG法) 不完全コレスキー分解による前処理CG法アルゴリズムの基本形を下図に 示す． 不完全コレスキー分解は元の行列要素の位置のみを計算する． ０要素に対するフィルインは計算しない． 計算は行列ベクトル積・ベクトルスカラ積・ベクトルの和・内積・除算で構成 される． 1.ICCG法のアルゴリズム 赤線で示した前処理部分は前進代入・後退代入で計算される． 2024年4月 ステップ1: αk = (rik •(LL T )−1 rik )/(Apik • p ik ) ステップ2: xi = xi +α pi ステップ3: ri k+ 1 = ri k − αk Apik ステップ4: β k = (rik+1 •(LLT ) −1 rik+1 ) /(ri k •(LLT ) −1 rik ) ステップ5: pi = (LL ) ri + β pi 計算科学技術特論B k+ 1 k+ 1 k k k T −1 k+1 22 k k

前進代入・後退代入処理 (L T )−1y = x T −1 (LL ) r = x (L T )−1L −1r = x L T は上三角行列である．Lと同様に 書き下す．前進代入処理と逆方向に 逐次に計算することによりxを求める −1 y = L r とする Ly = r と変形する L は下三角行列であるので 以下のように書ける． ことができる． (1)これを後退代入処理という L11y1 = r1 L21y1 + L22 y2 = r2 y1 = y2 = ⋯⋯ この式は右のように逐次に 解くことができる． (1)これを前進代入処理という 2024年4月 計算科学技術特論B r1 L11 1 (r − L21y1) L22 2 1 y3 = (r3 − L31y1 − L32 y2) L33 L31y1 + L32 y2 + L33 y3 = r3 L ⋯⋯ 23 後退代入の例 X 上三角行列 2 3 後退代入 7 9 10 後退代入 不完全コレスキー分解の前処理は本質的に リカレンスが発生する！ 2024年4月 計算科学技術特論B LT 24

不完全コレスキー分解のリカレンスの解消 不完全コレスキー分解の計算順序は前頁のようにリカレンスのために並列化できない． またガウス・ザイデル法と同様な参照関係を持つ(下図は前進代入の例)． そこで前述の(b)のハイバープレーン法や(c)のRED-BLACK法を用いループ内のリカレ ンス(依存関係)をなくす手法が使用できる． リカレンス(依存関係)をなくす事で前回講義のようにハードウェアの並列計算機構を使 い高速に計算可能となる．またコア間の並列化そのものにも使用可能． (a)オリジナルスイープ 参照関係 計算順序(スイープ順序) 2024年4月 計算科学技術特論B (b)ハイパープレーンスイープ 斜めのハイパープレーン内は リカレンスがない． ハイバープレーン間のみにリ カレンスが発生． (c)RED-BLACKスイープ 黄と青の２つのループに再構 成する． それぞれのループ内はリカレ ンスがない． 25 アプリケーションの 性能最適化 (高並列性能最適化) 2024年4月 計算科学技術特論B 26

計算科学 問題設定 計算機科学 理論に忠実な分かり やすいコーディング 理論構築 定式化 プログラム 離散化 現代のプログラミング (スパコンを使う場合) • アプリケーションの超 並列性を引き出す プログラム作成 • コンパイル デバッグ プロダクションRUN グラスィックス処理 言語 プログラム コンパイラ 結果の解釈 論文作成 評価 2024年4月 計算科学技術特論B 開発・実行環境 離散化 計算機科学 理論に忠実な分かり やすいコーディング 10万行以上 のプログラム • コンパイル デバッグ プロダクションRUN 評価 2024年4月 計算科学技術特論B プロセッサの単体性能を引き 出す 言語 プログラム 結果の解釈 論文作成 現代のプログラミング (スパコンを使う場合) • アプリケーションの超 並列性を引き出す プログラム作成 グラスィックス処理 システム運用 27 理論構築 定式化 ハードウェア 高並列化・ 高性能コーディング 計算科学 問題設定 プロセッサの単体性能を引き 出す コンパイラ 開発・実行環境 高並列化・ 高性能コーディング 28 ハードウェア システム運用

性能最適化技術の概要 高並列 単体性能 ソースコードの調査 現状認識 測定/評価法 計算・通信カーネルの決定 問題点把握 2024年4月 計算科学技術特論B 問題点の評価法 問題点の評価法 29 ソースコードの調査・測定/評価法・ 計算・通信カーネルの決定 高並列 単体性能 ソースコードの調査 現状認識 測定/評価法 計算・通信カーネルの決定 問題点把握 2024年4月 計算科学技術特論B 問題点の評価法 30 問題点の評価法

ソースコードの調査 測定/評価法 mainA subB subC subD DOループ subE 処理1 (1)コードの構造を分析し物理に沿った処理ブロック(計算/通信)に分割 (2)コードの実行時間の実測 (3)プログラムソースコードの調査 (4)処理ブロックの物理的処理内容を把握 (5)計算ブロック毎の計算特性把握(非並列/完全並列/部分並列、Nに比例/ N**2に比例等) (6)通信ブロック毎の通信特性把握 (グローバル通信・隣接通信、隣接面に比 例&隣接通信/体積に比例、等) 実行時間 物理的 演算・通信 演算・通信 ・スケーラ 処理内容 特性 見積り ビリティ ブロック1 (計算) 部分並列 Nに比例 処理1.1 DOループ 処理2 ブロック2 (計算) 処理2.1 処理3 (通信) 通信1 subF 2024年4月 完全並列 N**3に比例 ○ 隣接通信 隣接面に比例 ○ ブロック3 subG 計算・通信 カーネルの 決定 カー ネル (1)計算・通信ブロックについて物理的処理内容・コーディングの評価を行い同種の計算・ 通信ブロックを評価し異なる種類の計算・通信ブロックをカーネルの候補として洗い出す (2)並列特性分析の結果から得た問題規模に対する依存性の情報を元にターゲット問題 実行時に、また高並列実行時にカーネルとなる計算・通信カーネルを洗い出す 計算科学技術特論B 31 測定/評価法 ブロック毎の実行時間とスケーラビリティ評価(例) →従来の評価はサブルーチン毎・関数毎等の評価が多い →サブルーチン・関数は色々な場所で呼ばれるため正しい評価ができない LM+MSD法における測定区間毎のスケーラビリティ 10.0 全体 evolve_WFs_in_subspace 処理ブロック1 スケーラビリティ 処理ブロック２ m_ES_Vnonlocal_W 処理ブロック３ 処理ブロック４ modified_gram_schmidt 処理ブロック５ 処理ブロック６ m_XC_cal_potential betar_dot_WFs m_FFT_WF 1.0 128 2024年4月 計算科学技術特論B 256 並列数 32 512 1024

高並列における問題点の評価法 高並列 単体性能 ソースコードの調査 現状認識 測定/評価法 計算・通信カーネルの決定 問題点把握 2024年4月 問題点の評価法 計算科学技術特論B 評価法 問題点の評価法 33 ストロングスケールによる評価 １ノード 全体の問題規模を一定にして測定・評 価する方法． ここでは4×4×10が全体の問題規模． ２ノード ２ノード分割では１ノードあたりの問題規 模は4×4×5となる． 20ノード分割では１ノードあたりの問題 規模は2×2×2となる． 問題を１種類作れば良いので測定は楽 である． 20ノード 並列時の挙動は見えにくい． 2024年4月 計算科学技術特論B 34

ウィークスケーリングによる評価 評価法 １ノード ２ノード 1プロセッサあたりの問題規模を一定にして 測定・評価する方法． ここでは4×4×5が1プロセッサあたりの問 題規模． ２ノード分割では全体の問題規模は 4×4×10となる． 10ノード分割では全体の問題規模は 8×4×25となる． 問題を複数作る必要があり測定は煩雑で ある． ただし並列時の挙動が見え易い． 10ノード 2024年4月 計算科学技術特論B 35 35 大規模並列のウィークスケーリング評価 並列計算できない部分(1秒) 1+0.99=1.99秒 ほぼ50倍 並列計算できる部分(99秒) 1+0.099=1.099秒 ほぼ91倍 そもそも並列化とは？(3) 必要な並列度を確保する！ 1000 100 プロセッサ プロセッサ 2024年4月 計算科学技術特論B 並列計算できない部分(非並 列部)を限りなく小さくする！ 36

うとしても、アプリケーションの並列化の限界から、数千の並列度しか使えない、 している。しかし、このハードウェアの持つ 8 万並列以上の並列度を使お 万並列以上の並列度を使お としている。しかし、このハードウェアの持つ 8 万並列以上の並列度を使お 。しかし、このハードウェアの持つ 8 万並列以上の並列度を使お 場合が出てくる。これが、アプリケーションとハードウェアの並列度のミスマッチで プリケーションの並列化の限界から、数千の並列度しか使えない、という 列度しか使えない、という ションの並列化の限界から、数千の並列度しか使えない、という アプリケーションの並列化の限界から、数千の並列度しか使えない、という アプリケーションの並列度の限界に近づいてくると、演算時間が著しく小さくなる 並列度のミスマッチである。 。これが、アプリケーションとハードウェアの並列度のミスマッチである。 、アプリケーションとハードウェアの並列度のミスマッチである。 る。これが、アプリケーションとハードウェアの並列度のミスマッチである。 大規模並列のウィークスケーリング評価 し通信時間の割合が増大する状態となり、並列効率の悪化を招く。 が著しく小さくなるのに対 度の限界に近づいてくると、演算時間が著しく小さくなるのに対 ンの並列度の限界に近づいてくると、演算時間が著しく小さくなるのに対 ョンの並列度の限界に近づいてくると、演算時間が著しく小さくなるのに対 2 点目は、非並列部の残存である。本節の冒頭に述べたように、演算部分に非並 招く。 非並列部の評価 する状態となり、並列効率の悪化を招く。 合が増大する状態となり、並列効率の悪化を招く。 割合が増大する状態となり、並列効率の悪化を招く。 残っているとアムダールの法則により、並列性能の悪化を招く。ここで、あるアプ に、演算部分に非並列部が 残存である。本節の冒頭に述べたように、演算部分に非並列部が 並列部の残存である。本節の冒頭に述べたように、演算部分に非並列部が 非並列部の残存である。本節の冒頭に述べたように、演算部分に非並列部が ションの逐次実行時の実行時間を Ts であるとし、そのアプリケーションの並列化 。ここで、あるアプリケー の法則により、並列性能の悪化を招く。ここで、あるアプリケー ムダールの法則により、並列性能の悪化を招く。ここで、あるアプリケー アムダールの法則により、並列性能の悪化を招く。ここで、あるアプリケー であるとすると、アプリケーションの非並列化率は 1 − α となる。このアプリケー ケーションの並列化率を α 実行時間を T であるとし、そのアプリケーションの並列化率を α s 行時の実行時間を TsTであるとし、そのアプリケーションの並列化率を 実行時の実行時間を α α)) で s であるとし、そのアプリケーションの並列化率を を並列度 n で実行すると、 n 並列での実行時間 T は、 T = Ts ( αn + (1α− n n る。このアプリケーション ケーションの非並列化率は 1 − α となる。このアプリケーション 、アプリケーションの非並列化率は 1− αα となる。このアプリケーション と、アプリケーションの非並列化率は 1− となる。このアプリケーション αる。仮に n = 10000 の時に並列化効率α50%を達成するための α を求めると、並列 T ( + (1 − α)) で表され s nn 並列での実行時間 Tn は、Tn = Ts ( + (1 − α)) と、 n nTn==TsT 実行すると、 nを要求されることになる。非並列部の残存を発見しやすいのは、前述し 並列での実行時間 +(1(1−−α)) α))で表され で表され 行すると、 n 並列での実行時間 TnTは、 T (sαn( αn+で表され n は、 は 99.99% α を求めると、並列化率 α 時に並列化効率 50%を達成するための α を求めると、並列化率 α 10000 の時に並列化効率 50% を達成するための を求めると、並列化率αα 000 の時に並列化効率 50% を達成するための αα を求めると、並列化率 ウィークスケーリング測定による、演算部の実行時間の増大である。 やすいのは、前述した通り ことになる。非並列部の残存を発見しやすいのは、前述した通り α 3 点目は、大域通信における大きな通信サイズ、通信回数の発生である。通信時 求されることになる。非並列部の残存を発見しやすいのは、前述した通り されることになる。非並列部の残存を発見しやすいのは、前述した通り T = T + Ts (1− α ) 並列度を上げると，ここが である。 n s 定による、演算部の実行時間の増大である。 どんどん大きくなる n に大域通信は、並列性能に大きな影響を与える。n ノード間で m[Byte] の allredu ーリング測定による、演算部の実行時間の増大である。 リング測定による、演算部の実行時間の増大である。 発生である。通信時間、特 おける大きな通信サイズ、通信回数の発生である。通信時間、特 α を実施する例を考える。2α分木のアルゴリズムにより allreduce 通信を実施し通信 大域通信における大きな通信サイズ、通信回数の発生である。通信時間、特 域通信における大きな通信サイズ、通信回数の発生である。通信時間、特 T = 2T + 2T (1− α ) = T + 2T (1− α) で m[Byte] の allreduce 通信 2n s s s s 能に大きな影響を与える。 n ノード間で m[Byte] の allreduce 通信 P [B/sec] と仮定すると、全ノードが通信を終わり m[Byte] の総和量を取得するた 2n n ノード間で nm[Byte]ののallreduce t 、並列性能に大きな影響を与える。 n ノード間で allreduce通信 通信 並列性能に大きな影響を与える。 m[Byte] ce 通信を実施し通信性能が m×log 2nallreduce 通信を実施し通信性能が 。2 分木のアルゴリズムにより 信時間 T は、T = で計算される。大域通信と隣接通信と比較するために、 Pt を考える。 2 分木のアルゴリズムにより allreduce 通信を実施し通信性能が 考える。 2 分木のアルゴリズムにより allreduce 通信を実施し通信性能が 総和量を取得するための通 、全ノードが通信を終わり m[Byte] の総和量を取得するための通 ドが隣のランクに m[Byte] の隣接通信を実施する例を考える。通信性能を上記の条 定すると、全ノードが通信を終わり m[Byte] の総和量を取得するための通 m 2n 信と比較するために、 n ノー すると、全ノードが通信を終わり m[Byte] の総和量を取得するための通 で計算される。大域通信と隣接通信と比較するために、 n ノー わせると、全ノードが m[Byte] の通信が終了するまでの通信時間 T は、T = Pt で m×log 2n m×log 2n T = で計算される。大域通信と隣接通信と比較するために、 ノー 通信性能を上記の条件と合 =e] n nノー の隣接通信を実施する例を考える。通信性能を上記の条件と合 Pt Pt で計算される。大域通信と隣接通信と比較するために、 れる。大域通信と隣接通信を比較すると、 log 2n の係数分だけ大域通信時間の方が m 時間 Tm[Byte] は、 Tの隣接通信を実施する例を考える。通信性能を上記の条件と合 = で計算さ クに の隣接通信を実施する例を考える。通信性能を上記の条件と合 に m[Byte] P t m[Byte] の通信が終了するまでの通信時間 T は、T = Pmt で計算さ ことがわかる。大域通信は最低限必要な量だけを通信するようにすべきである。 n 大規模並列のウィークスケーリング評価 m m 大域通信時間の方が大きい ノードが m[Byte] の通信が終了するまでの通信時間 は、 で計算さ ードが m[Byte] の通信が終了するまでの通信時間 T Tは、 T T==PtPで計算さ 信を比較すると、 log 2n の係数分だけ大域通信時間の方が大きい t の計算体系で各ノードが m 要素のベクトル量を保持している場合に、全ノードで うにすべきである。 n ノード loglog n 現状使用可能な実行環境を使用し100程度/1000程度/数千程度と段階を追っ 信と隣接通信を比較すると、 の係数分だけ大域通信時間の方が大きい と隣接通信を比較すると、 2n2n の係数分だけ大域通信時間の方が大きい は最低限必要な量だけを通信するようにすべきである。 n ノード 要素の総和を計算する例を考える。時々見受けられる例であるが、全ノードで m要 てできるだけ高い並列度で並列性能を確認する(ウィークスケーリング測定). 場合に、全ノードで n×m 。大域通信は最低限必要な量だけを通信するようにすべきである。 ノード クトル量を allreduce しその後、各ノード内で m 要素の和を計算している場合があ が m 要素のベクトル量を保持している場合に、全ノードで n × m難しい場合 n ウィークスケーリングが難しいものもあるが出来るだけ測定したい. 大域通信は最低限必要な量だけを通信するようにすべきである。 n nノード が、全ノードで m 要素のベ は，演算時間をモデル化して実測とモデルとの一致を評価する. 各ノードが 要素のベクトル量を保持している場合に、全ノードで を考える。時々見受けられる例であるが、全ノードで m 要素のベ ノードが mm 要素のベクトル量を保持している場合に、全ノードで n n××mm 算している場合がある。こ 計算する例を考える。時々見受けられる例であるが、全ノードで 要素のベ 計算のブロック毎に実施 その後、各ノード内で m 要素の和を計算している場合がある。こ 29 算する例を考える。時々見受けられる例であるが、全ノードで mm要素のベ 非並列部の評価 reduce しその後、各ノード内で 要素の和を計算している場合がある。こ duce しその後、各ノード内で mm 要素の和を計算している場合がある。こ 2024年4月 計算科学技術特論B 37 演算時間が増大しないか確認．演算時 間が増大すれば非並列部が残っている 29演算時間 2929 実行 時間 並列数 2024年4月 計算科学技術特論B 38

並列化に使用する通信 計算時間のバラつき =ロードインバランス 計算時間 通信時間 通信時間 プロセッサ１ 隣接通信 プロセッサ２ 隣接通信 プロセッサ３ 隣接通信 大域通信 例えば総和計算 MPI_allreduce プロセッサ４ 2024年4月 計算科学技術特論B 39 MPIの概要 2024年4月 計算科学技術特論B 40

超高並列を目指した場合の評価点 −ブロック毎に以下を評価するn n n n 非並列部が残っていないか？残っている場合に問題ないか？ 隣接通信時間が超高並列時にどれくらいの割合を占めるか？ 大域通信時間が超高並列時にどれくらい増大するか？ ロードインバランスが超高並列時に悪化しないか？ これらの評価が重要 そのために 2024年4月 計算科学技術特論B 41 ストロングスケールとウィークスケーリング 測定を使う評価 n ターゲット問題を決める→どれくらいの問題を解きたいか？そのためにどれく らいのノード数が必要か？ n １プロセッサの問題規模がターゲット問題と同程度でできればウィークスケー リングで実行時間・ロードインバランス・隣接通信時間・大域通信時間を測定・ 評価する n 難しければ100ノード並列程度までストロングスケールで測定しても良い n その程度の並列度でスケールしていなければ根本的な問題があるので解決する n 問題なければ並列度を上げてウィークスケーリングで大規模並列の挙動を測 定する ストロング スケール 実行 時間 ウィークスケール 1プロセッサの問題の大きさをターゲット問題の大きさに合わせる 並列数 2024年4月 計算科学技術特論B 42

大規模並列のウィークスケーリング評価 n 現状使用可能な実行環境を使用し100程度/1000程度/数千程度と段階を追っ てできるだけ高い並列度で並列性能を確認する(ウィークスケーリング測定). n ウィークスケーリングが難しいものもあるが出来るだけ測定したい. 難しい場合 は，演算時間をモデル化して実測とモデルとの一致を評価する. 計算のブロック毎に実施 非並列部の評価 演算時間が増大しないか確認．演算時 間が増大すれば非並列部が残っている 演算時間 大域通信時間の割合は一般に増大する. 通信量・通信回数を確認する事で増大の割合を見極める. 実行 時間 大域通信時間 隣接通信時間 通信時間の評価 隣接時間の割合と増大しないか確認. 通信の処理方法が悪いと隣接通信時間が増大する. 並列数 2024年4月 計算科学技術特論B 43 Weak Scaling評価例 処理ブロックA 処理ブロックB B A 処理ブロックC 処理ブロックD D C 2024年4月 計算科学技術特論B 44

1000Proc. Weak Scaling測定例 (横軸をランク番号として評価) 隣接通信 大域通信 通信時間の評価 in grid points are kept within a part and do not 演算 affect the communications within the other part. The number of parts is taken up to the number of parallelisms in orbitals. We used this mapping in the experiments. Data-1 on 12288, 24576, 26864 and 73728 cores (cases 1, 2, 3 and 4, respectively). The number of parallel tasks is the product of the number of parallel tasks in space and that in orbitals. The number of parallel tasks in space was fixed to 1,536, and the number of parallel tasks in orbitals was varied from 1, 2, 3 and 6. Therefore, the number of プロセス番号 parallel tasks amounted to 1536, 3072, 4608, or 9216. 5. PERFORMANCE ANALYSIS 8000Proc. The scalability and efficiency of the modified code which is parallelized in terms of grid points and orbitals as described above were evaluated with a simulation model named “Data-1,” depicting a SiNW with 19,848 atoms. The number of grid 隣接通信 points and orbitals were 320x320x120 and 40,992, respectively. We used an MPI library tuned to the Tofu network to obtain higher 大域通信 performance. We also made an appropriate task mapping to the Tofu network. Figure 6 plots the computation time and communication time for GS, CG, MatE in SD (MatE/SD) and RotV in SD (RotV/SD) . The horizontal axis in each figure stands for the number of cores. The theoretical computation time for cases 2, 3 and 4 is also ploted where the computation time of case1 is divided by 2, 3 and 6, respectively. They show the scalability in the parallelization for the orbitals. The measured computation time becames closer to the theoretical computation time as the number of parallel tasks in orbitals increased, and GS, CG, MatE/SD and RotV/SD were found to be well-parallelized in orbitals. ロードインバランス の評価 5.1 Performance of DGEMM The principal kernel of the modified RSDFT code is a DGEMM routine of the BLAS. We optimized the DGEMM routine so as to utilize the special features such as SIMD instructions and the sector cache and hardware barrier synchronization capability of the SPARC64TM VIIIfx effectively. The sustained performance of DGEMM on a compute node is 123.7 giga-flops, or 96.6% of the peak performance. In particular, we found that the computation 2024年4月 計算科学技術特論B time as a result of keeping the block data on the L1 cache manually decreased by 12% compared with the computation time for the usual data replacement operations of the L1 cache. This DGEMM tuned for the K computer was also used for the LINPACK benchmark program. 5.2 演算 45 The communication time with respect to the parallel tasks in space decreased in proportion to the number of parallel tasks in orbitals, because the number of calls of MPI functions for global and adjacent communications decreased in proportion to the reciprocal of the number of parallel tasks in orbitals. On the other hand, the global communication time for the parallel tasks in orbitals was supposed to increase as the number of parallel tasks in orbitals increased. The number of MPI processes requiring communications for the parallel tasks in orbitals, however, was actually restricted to a relatively small number of compute nodes, and therefore, the wall clock time for global communications of the parallel tasks in orbitals was small. This means we succeeded in decreasing time for global communication by the combination ストロングスケーリングだけどスケールを Scalability モデル化した例 We measured the computation time for the SCF iterations with (a) (b) 400.0 160.0 theoretical computation theoretical computation computation adjacent/space global/space global/orbital computation global/space 120.0 global/orbital Time per CG (sec.) Time per GS (sec.) 300.0 wait/orbital 200.0 80.0 40.0 100.0 0.0 0.0 0 20,000 40,000 60,000 Number of cores (c) 0 80,000 60000 80000 300.0 150.0 theoretical computation theoretical computation computation computation adjacent/space adjacent/space global/space Time per RotV/SD(sec.) Time per MatE/SD (sec.) 40000 Number of cores (d) 200.0 global/orbital 100.0 50.0 global/space global/orbital 200.0 100.0 0.0 0.0 0 2024年4月 20000 20,000 計算科学技術特論B 40,000 Number of cores 60,000 0 80,000 46 20,000 40,000 Number of cores 60,000 Figure 6. Computation and communication time of (a) GS, (b) CG, (c) MatE/SD and (d) RotV/SD for different numbers of cores. 80,000

並列性能上のボトルネック 今まで示した評価を実施することにより処理ブロック毎に並列性能上の 問題がある事が発見される． それらを分析するとだいたい以下の６点に分類されると考える． 2024年4月 1 アプリケーションとハードウェアの並列度のミスマッチ (アプリケーションの並列度不足) 2 非並列部の残存 3 大域通信における大きな通信サイズ、通信回数の発生 4 フルノードにおける大域通信の発生 5 隣接通信における大きな通信サイズ、通信回数の発生 6 ロードインバランスの発生 計算科学技術特論B 47 アプリケーションとハードウェアのミスマッチ 波数展開を使う第一原理計算の場合を考 える． 電子のエネルギーバンド並列のみ実装され ていると仮定 原子一個に対し最外核の電子2個と考える とエネルギーバンドの数は，おおよそ原子 数：N×２となる． N=10000ならエネルギーバンド数は 20000程度である． この場合，最大でも20000並列までしか到 達しない． 通信の増大や計算の粒度を考えると数百 並列で限界となる． 数万のノードを装備するスパコンとアプリ ケーションの並列度がミスマッチしている 波数等その他の分割を組み合わせた並列 化が必要となる． 詳細は3回目の講義で実例を示す． 2024年4月 2022年4月 計算科学技術特論B 48

! ! ! 非並列部の残存 BLAS Level3 BLAS Level3 ( 2) ( ( 第一原理アプリケーションの処理ブロック:ノンローカル項の計算部分を例に示す． ( 2) 低並列でストロングスケールで測定． すでにこの並列度でスケールしていない．原因は次ページに示す． 14 [s] 1 12 1 10 1 8 6 4 2 0 16 2024年4月 計算科学技術特論B 49 非並列部の残存 第一原理計算のアプリでエネルギーバンド並列を使用している場合の非並列部が残っ HPCS2012 HPCS2012 ている例． subroutine m_es_vnonlocal_w(ik,iksnl,ispin,switch_of_eko_part) | +-call tstatc0_begin | | loop_ntyp: do it = 1, ntyp | | | loop_natm : do ia = 1, natm ---------原子数のループ | | | +-call calc_phase | | | | T-do lmt2 = 1, ilmt(it) | | | | +-call vnonlocal_w_part_sum_over_lmt1 | | | | | +-call add_vnlph_l_without_eko_part | | | | | | subroutine add_vnlph_l_without_eko_part() | | | | | | | T-if(kimg == 1) then | | | | | | | | T-do ib = 1, np_e ---------エネルギーバンド並列部 | | | | | | | | | T-do i = 1, iba(ik) | | | | | | | | | V-end do | | | | | | | | V-end do | | | | | | | +-else | | | | | | | | T-do ib = 1, np_e ---------エネルギーバンド並列部 | | | | | | | | | T-do i = 1, iba(ik) | | | | | | | | | V-end do | | | | | | | | V-enddo | | | | | | | V-end if | | | | | | end subroutine add_vnlph_l_without_eko_part | | | | V-end do | | | V-end do loop_natm | | V-end do loop_ntyp end subroutine m_es_vnonlocal_w 2024年4月 計算科学技術特論B 50

大域通信における大きな通信サイズ・通信回数 各ノードでベクトルデータを総和しスカラー値を求め，スカラー値をノード 間でallreduceすれば良いのに，全ノードで巨大なベクトルをallreduceして いませんか？ この場合は，よけいなバンド幅を消費していることになる． (1)総和 (2)総和 計算 計算 P1 P1 P2 P2 ・・・・・・・ ・・・・・・・ Pi Pi (1)ALLREDUCE ・・・・・・・ ・・・・・・・ Pn 2024年4月 (2)ALLREDUCE Pn 計算科学技術特論B 51 大域通信における大きな通信サイズ・通信回数 Efficiency of Mapping to the Tofu network - Gram-Schmidt Optimal Mapping ! avoid some conflict among sub-communicators 通信はレイテンシーが問題になる場合とバンド selected Tofu specific communication algorithm 幅が問題になる場合がある． • • • • 120.0## 100.0## Time (sec.) 80.0## MPI_Bcast 60.0## 40.0## MPI_Allreduce 20.0## 0.0## Communication Wati/Orbital 1D 12.929 Communication Grobal/Orbital 17.118 Communication Global/Space 33.565 Computation 32.050 2024年4月 計算科学技術特論B 各ノードでベクトルデータを総和しスカラー値を SiNW 19,848 atoms 求め，スカラー値をノード間でallreduceすれば Grids : 320x320x120 良いのに，全ノードで巨大なベクトルを Orbitals: 41,472 allreduceしていませんか？ Total process number: 12,288 この場合は，よけいなバンド幅を消費している " Grid space: 2,048(32x32x2) " Orbital: 6 ことになる． • Torus Mapping :32x32x12 最適化された大域通信のライブラリを使ってい ますか？ 最適なライブラリと最適化されていないライブ ラリでは数倍性能差が発生する場合がある． 小さなデータを頻繁に大域通信していません か？ Optimal 3D まとめることが出来る通信はまとめて回数を減 8.295 らす事でレイテンシーネックを解消することがで 3.726 きます． 3.957 大域通信ライブラリで可能な通信と同様な事 32.111 を隣接通信で自作しているような場合がある． 性能の良い大域通信ライブラリの使用を推奨 します． 52

い並列効果を得るには、オリジナルコードを使用する限り、ノードあたりの原 と大きくとる必要がある。したがって区間 5 でも、区間 2 と同様のハードウェ のアプリの並列度のミスマッチが発生していると言える。Fig. 4.25 の通信処理 節の Fig. 4.22 に示したトランスバース転送である。スケーラビリティの悪化の は、Fig. 4.22 に示したトランスバース転送による大域通信時間の増大である。 フルノードでAlltoallのような通信を実施す 間自体も良いスケーラビリティを示しているわけではないため、演算自体にも る事は性能を著しく劣化させる． 化の原因が生じているものと考えられる。 第一原理アプリケーションのそのような通 フルノードにおける大域通信の発生 信を含む処理ブロックの例を左図に示す． すでに低並列度でストロングスケールして いない． 並列軸の追加によりFFTのAlltoall通信を フルノードでなく一部のノードに限ることが 出来る場合がある． 詳細は3回目の講義で説明する． CG法の内積計算のようにフルノードで大 域通信を実施する必要がある計算は減ら すことができない． その場合は，高速なハードウェアアシスト 機能を使用する等の対処をとる． 計算科学技術特論B 53 Fig. 4.25: PHASE 区間 5 のスケーラピリティ 2024年4月 隣接通信における大きな通信サイズ・通信回数 む処理 隣接通信がネックになる場合，アルゴリズムを見直す事で隣接 FFT を含む処理となっており、基本的にエネルギーバンド並列で並列化され 通信をやめてしまうことが出来る場合がある． 例えばCG法前処理のローカライズ化等の場合． 3.3 Method for Parallelization HfSiO2、384 原子アモルファス系を入力データとして、 16 から 128 並列の低 トロングスケールで測定した。この計算体系のエネルギーバンド数は 1280 で 中途半端な並列化が通信性能の悪化を招く場合がある． 1)Forward/Backward Substitution これらの詳細は3回目の講義で説明する． 1280 のエネルギーバンド数を 16 から 128 並列で並列化していることとなる。 ->Localized ILU Method リティのデータを Fig. 4.26 に示す。この図の中で、目立っている２つの系列 (b)ハイパープレーンスイープ うな処理を実行している。 PE1 plane2 plane1 PE (1) plane3 述べた行列行列積にできる処理 PE2 (2) cal cal cal com cal com cal PE1 (1) cal PE2 (2) PE3 (3) PE4 (4) (3) com com cal cal cal バンド数分のPE3FFT を含む計算処理 PE4 (4) cal com cal com PE cal plane1 plane2 cal cal plane3 ・・・・・ cal cal com cal com cal cal com cal com cal cal com cal com cal PE5 (5) の中で系列 PE5 2 については、行列行列積にできる処理で述べたように、 64 から (5) com com cal cal com cal com cal cal cal (6) 全く性能向上していない。系列 の FFT PE6 を含む処理ついてはバンド計算の中で com com PE6 (6) cal com 5 cal cal cal com cal cal TIME TIME T を実行しており、バンドによる並列は実行されているのでよくスケールして cal : calculation time cal : calculation time com : latency time of communication com : latency time of communication 見える。ただし系列 2・5 を含む処理時間全体を見ると、ノードあたり 80 程度 を分担させたとしても、系列 2 の処理が原因となって並列のスケーラビリティ 3)Process of CG Method 計算科学技術特論B ->There is no ploblem. 2024年4月 54

ロードインバランスの発生 ロードインバランスの発生は色々な原因がある． 第一は粒子の計算で時間発展により粒子のバランスが変わるような場合． . . . . . 動的負荷分散を取る方法等色々な手法が研究されている． NICAM gl07 rl00 gl14 rl07 その他システムの外乱が原因となる場合もある． phystep,dynstep 詳細は３回目以降に講義する． 以下の例は処理ブロック毎にスケーラビリティを測定した例． . 2024年4月 gl7 計算科学技術特論B gl14 -1 NICAM phystep 81920 rd_driver gl13 20 55 rd_driver 簡単な実例 高並列 単体性能 ソースコードの調査 現状認識 . 測定/評価法 -2 NICAM dynstep dynstep gl14 81920 計算・通信カーネルの決定 問題点把握 2024年4月 計算科学技術特論B 問題点の評価法 56 問題点の評価法

ソースコードの調査 測定/評価法 mainA subB subC subD DOループ subE 処理1 (1)コードの構造を分析し物理に沿った処理ブロック(計算/通信)に分割 (2)コードの実行時間の実測 (3)プログラムソースコードの調査 (4)処理ブロックの物理的処理内容を把握 (5)計算ブロック毎の計算特性把握(非並列/完全並列/部分並列、Nに比例/ N**2に比例等) (6)通信ブロック毎の通信特性把握 (グローバル通信・隣接通信、隣接面に比 例&隣接通信/体積に比例、等) 実行時間 物理的 演算・通信 演算・通信 ・スケーラ 処理内容 特性 見積り ビリティ ブロック1 (計算) 部分並列 Nに比例 処理1.1 DOループ 処理2 ブロック2 MG_old 20240408T172529 mg.f (564 nodes) (計算) 処理2.1 call mpi_init(ierr)[mg_mpi:mg.f:92-92] 処理3 call mpi_comm_rank(mpi_comm_world, me, ierr)[mg_mpi:mg.f:93-93] (通信) call mpi_comm_size(mpi_comm_world, nprocs, ierr)[mg_mpi:mg.f:94-94] 通信1 program mg_mpi[mg_mpi:mg.f:49-396] カー ネル 完全並列 N**3に比例 ○ 隣接通信 隣接面に比例 ○ if (nprocs_compiled .gt. maxprocs) then[mg_mpi:mg.f:97-103] subF do i = 1, t_last[mg_mpi:mg.f:112-114] ブロック3 call mpi_barrier(MPI_COMM_WORLD, ierr)[mg_mpi:mg.f:116-116] subG call timer_start(T_init)[mg_mpi:mg.f:118-118] if( me .eq. root )then[mg_mpi:mg.f:125-155] call mpi_bcast(lt, 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:157-157] call mpi_bcast(nit, (1)計算・通信ブロックについて物理的処理内容・コーディングの評価を行い同種の計算・ 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:158-158] 計算・通信 通信ブロックを評価し異なる種類の計算・通信ブロックをカーネルの候補として洗い出す call mpi_bcast(nx(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:159-159] カーネルの (2)並列特性分析の結果から得た問題規模に対する依存性の情報を元にターゲット問題 call mpi_bcast(ny(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:160-160] 決定call mpi_bcast(nz(lt), 1,実行時に、また高並列実行時にカーネルとなる計算・通信カーネルを洗い出す MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:161-161] call mpi_bcast(debug_vec(0), 8, MPI_INTEGER, 0,[mg_mpi:mg.f:162-163] 2024年4月 call 計算科学技術特論B mpi_bcast(timeron, 1, MPI_LOGICAL, 0, mpi_comm_world, ierr)[mg_mpi57 :mg.f:164-164] call setup(n1,n2,n3,k)[mg_mpi:mg.f:224-224] subroutine setup(n1,n2,n3,k)[setup:mg.f:401-565] call zero3(u,n1,n2,n3)[mg_mpi:mg.f:225-225] CCA/EBT call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:226-226] MG_old 検索 (Code Comprehension Assistance for EvidenceBased[ performance Tuning) :mg.f:249-249] call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:228-228] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:248-248] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) mg_mpi 20240408T172529 call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:256-256] mg.f (564 nodes) call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:257-257] program mg_mpi[mg_mpi:mg.f:49-396] call setup(n1,n2,n3,k)[mg_mpi:mg.f:258-258] call mpi_init(ierr)[mg_mpi:mg.f:92-92] call zero3(u,n1,n2,n3)[mg_mpi:mg.f:259-259] call mpi_comm_rank(mpi_comm_world, me, ierr)[mg_mpi:mg.f:93-93] call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:260-260] call mpi_comm_size(mpi_comm_world, nprocs, ierr)[mg_mpi:mg.f:94-94] call timer_stop(T_init)[mg_mpi:mg.f:262-262] if (nprocs_compiled .gt. maxprocs) then[mg_mpi:mg.f:97-103] do i = 1, t_last[mg_mpi:mg.f:269-271] do i = 1, t_last[mg_mpi:mg.f:112-114] call mpi_barrier(mpi_comm_world,ierr)[mg_mpi:mg.f:272-272] call mpi_barrier(MPI_COMM_WORLD, ierr)[mg_mpi:mg.f:116-116] call timer_start(T_bench)[mg_mpi:mg.f:274-274] call timer_start(T_init)[mg_mpi:mg.f:118-118] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:276-276] if( me .eq. root )then[mg_mpi:mg.f:125-155] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:277-277] call mpi_bcast(lt, 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:157-157] do it=1,nit[mg_mpi:mg.f:281-288] call mpi_bcast(nit, 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:158-158] call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:286-286] call mpi_bcast(nx(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:159-159] subroutine mg3P(u,v,r,a,c,n1,n2,n3,k)[mg3p:mg.f:570-632] call mpi_bcast(ny(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:160-160] do k= lt, lb+1 , -1[mg3p:mg.f:594-598] call mpi_bcast(nz(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:161-161] call rprj3(r(ir(k)),m1(k),m2(k),m3(k),[mg.f:596-597] call mpi_bcast(debug_vec(0), 8, MPI_INTEGER, 0,[mg_mpi:mg.f:162-163] subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k )[rprj3:mg.f:777-869] call mpi_bcast(timeron, 1, MPI_LOGICAL, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:164-164] if (timeron) call timer_start(t_rprj3)[rprj3:mg.f:803-803] call setup(n1,n2,n3,k)[mg_mpi:mg.f:224-224] do j3=2,m3j-1[rprj3:mg.f:822-853]Br:0,St:13,FOp:24, ES2 AR2:5,DAR2:5,IAR2:0,B/F:1.67 subroutine setup(n1,n2,n3,k)[setup:mg.f:401-565] do j2=2,m2j-1[rprj3:mg.f:825-852] call zero3(u,n1,n2,n3)[mg_mpi:mg.f:225-225] do j1=2,m1j[rprj3:mg.f:829-836] call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:226-226] do j1=2,m1j-1[rprj3:mg.f:838-850] call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:228-228] if (timeron) call timer_stop(t_rprj3)[rprj3:mg.f:854-854] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:248-248] call comm3(s,m1j,m2j,m3j,j)[mg.f:858-858] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:249-249] if( debug_vec(0) .ge. 1 )then[rprj3:mg.f:860-862] call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:256-256] if-then-block[rprj3:mg.f:860-861] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:257-257] call rep_nrm(s,m1j,m2j,m3j,' rprj3',k-1)[mg.f:861-861] call setup(n1,n2,n3,k)[mg_mpi:mg.f:258-258] if( debug_vec(4) .ge. k )then[rprj3:mg.f:864-866] call zero3(u,n1,n2,n3)[mg_mpi:mg.f:259-259] call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:604-604] call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:260-260] call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)[mg.f:605-605] call timer_stop(T_init)[mg_mpi:mg.f:262-262] do k = lb+1, lt-1[mg3p:mg.f:607-623] 2024年4月 do 計算科学技術特論B 58 i = 1, t_last[mg_mpi:mg.f:269-271] call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:612-612] call mpi_barrier(mpi_comm_world,ierr)[mg_mpi:mg.f:272-272] call interp(u(ir(j)),m1(j),m2(j),m3(j),[mg.f:613-614] call timer_start(T_bench)[mg_mpi:mg.f:274-274] subroutine interp( z,mm1,mm2,mm3,u,n1,n2,n3,k )[interp:mg.f:874-1038] プログラム構造を俯瞰的にみて 理解を助けるツール． ・・・・・・ Fortran及びCが対象． プロシージャ呼び出し,ループ,分 岐に代表されるプログラムの構造 の可視化を行う． Loop Type プログラム構造解析支援ツール． 多数のFortranプログラムの解析 実績あり．

NPB MGを使用した解析例 NPB MG(マルチグリッド法により連立一次方程式 を解くベンチマークテストプログラム) Vサイクルマルチグリッド ファイングリッド コーズグリッド 2024年4月 計算科学技術特論B 59 NPB MGを使用した解析例 NPB MG(マルチグリッド法により連立一次方程式 を解くベンチマークテストプログラム) MGのアルゴリズム Au = f を解く Vサイクルマルチグリッド ファイングリッド Ak+1uk+1 = rk+1 制限補間 rk+1 = fk+1 − Ak+1uk+1 fk = Rk+1→k fk+1 fk−1 = Rk→k−1 fk uk+1 = Rk→k+1uk fk−2 = Rk−1→k−2 fk−1 Ak uk = rk 近似解求解 rk = fk − Ak uk 残差計算 f1 = R2→1 f2 uk = Rk−1→k uk−1 延長補間 コーズグリッド A1u1 = f1 2024年4月 計算科学技術特論B 60

call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:226-226] call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:228-228] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:248-248] 529 アプリケーションの構造の調査(概要) call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:249-249] call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:256-256] (564 nodes) call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:257-257] MG_old ogram mg_mpi[mg_mpi :mg.f:49-396] call setup(n1,n2,n3,k)[mg_mpi:mg.f:258-258] Kernel? Code < > Clear call zero3(u,n1,n2,n3)[mg_mpi:mg.f:259-259] call mpi_init(ierr)20240408T172529 [mg_mpi:mg.f:92-92] (1) 計測のスタートとストップに注目 call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:260-260] mg.f (564 nodes) call mpi_comm_rank(mpi_comm_world, me, ierr)[mg_mpi:mg.f:93-93] call timer_stop(T_init)[mg_mpi:mg.f:262-262] program mg_mpi[mg_mpi:mg.f:49-396] call mpi_comm_size(mpi_comm_world, nprocs, ierr)[mg_mpi:mg.f:94-94] do i = 1, t_last[mg_mpi :mg.f:269-271] ・・・・ call mpi_init(ierr)[mg_mpi:mg.f:92-92] call mpi_barrier(mpi_comm_world,ierr)[mg_mpi:mg.f:272-272] if (nprocs_compiled .gt. maxprocs) then[mg_mpi:mg.f:97-103] (2) 繰り返しループに注目 (3) mg3pサブルーチンに注目 call mpi_comm_rank(mpi_comm_world, me, ierr)[mg_mpi:mg.f:93-93] call timer_start(T_bench) [mg_mpi:mg.f:274-274] do i = 1, t_last[mg_mpi:mg.f:112-114] call mpi_comm_size(mpi_comm_world, nprocs, ierr)[mg_mpi:mg.f:94-94] call resid(u,v,r,n1,n2,n3,a,k) [mg_mpi:mg.f:276-276] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) [mg_mpi :mg.f:277-277] if (nprocs_compiled .gt.:mg.f:116-116] maxprocs) then[mg_mpi :mg.f:97-103] call mpi_barrier(MPI_COMM_WORLD, ierr)[mg_mpi do it=1,nit :mg.f:281-288] do i [=mg_mpi 1, t_last [mg_mpi:mg.f:112-114] call timer_start(T_init)[mg_mpi:mg.f:118-118] (4) マイナス方向のループでダウン サイクルと予測 call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:286-286] call mpi_barrier(MPI_COMM_WORLD, ierr)[mg_mpi:mg.f:116-116] if( me .eq. root )then[mg_mpi:mg.f:125-155] subroutine mg3P(u,v,r,a,c,n1,n2,n3,k)[mg3p:mg.f:570-632] call timer_start(T_init)[mg_mpi:mg.f:118-118] (5) 制限補間ルーチンと予測 do mpi_comm_world, k= lt, lb+1 , -1[mg3p:mg.f:594-598] call mpi_bcast(lt, 1, MPI_INTEGER, 0, ierr)[mg_mpi:mg.f:157-157] if( me .eq. root )then[mg_mpi:mg.f:125-155] call rprj3(r(ir(k)),m1(k),m2(k),m3(k),[mg.f:596-597] (6) ダウンサイクル後の近似解求解 と予測 call mpi_bcast(lt, 1, MPI_INTEGER, 0, [mpi_comm_world, ierr)[mg_mpi:mg.f:157-157] call mpi_bcast(nit, 1, MPI_INTEGER, 0, mpi_comm_world, ierr) mg_mpi:mg.f:158-158] call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:604-604] call mpi_bcast(nx(lt), 1, call mpi_bcast(nit, 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:158-158] MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi :mg.f:159-159] call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k) [mg.f:605-605] call mpi_bcast(nx(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:159-159] do k = lb+1, lt-1[mg3p:mg.f:607-623] call mpi_bcast(ny(lt), 1, MPI_INTEGER, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:160-160] (7) プラス方向のループでアップサ イクルと予測 call mpi_bcast(ny(lt), 1, MPI_INTEGER, 0, [mg.f:612-612] mpi_comm_world, ierr)[mg_mpi:mg.f:160-160] call zero3(u(ir(k)),m1(k),m2(k),m3(k)) call mpi_bcast(nz(lt), 1, MPI_INTEGER, 0, interp(u(ir(j)),m1(j),m2(j),m3(j), mpi_comm_world, ierr) [mpi_comm_world, mg_mpi:mg.f:161-161] call [mg.f:613-614] call mpi_bcast(nz(lt), 1, MPI_INTEGER, 0, ierr)[mg_mpi:mg.f:161-161] call resid(u(ir(k)),r(ir(k)),r(ir(k)),m1(k),m2(k),m3(k),a,k) mpi_bcast(debug_vec(0), 8, MPI_INTEGER, 0,[mg_mpi[mg.f:618-618] :mg.f:162-163] call mpi_bcast(debug_vec(0), 8, call MPI_INTEGER, 0,[mg_mpi :mg.f:162-163] call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)[mg.f:622-622] (8) 延長補間と予測 call mpi_bcast(timeron, 1, MPI_LOGICAL, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:164-164] call mpi_bcast(timeron, 1, MPI_LOGICAL, 0, mpi_comm_world, ierr)[mg_mpi:mg.f:164-164] call interp(u(ir(j)),m1(j),m2(j),m3(j),u,n1,n2,n3,k)[mg.f:627-627] call call setup(n1,n2,n3,k)[mg_mpi:mg.f:224-224] call resid(u,v,r,n1,n2,n3,a,k)[mg.f:628-628] setup(n1,n2,n3,k)[mg_mpi:mg.f:224-224] subroutine setup(n1,n2,n3,k)[setup:mg.f:401-565] call psinv(r,u,n1,n2,n3,c,k)[mg.f:629-629] subroutine setup(n1,n2,n3,k)call [setup :mg.f:401-565] :mg.f:225-225] zero3(u,n1,n2,n3)[[mg_mpi call resid(u,v,r,n1,n2,n3,a,k) mg_mpi:mg.f:287-287] call zero3(u,n1,n2,n3)[mg_mpi :mg.f:225-225] call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi :mg.f:226-226] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) [mg_mpi :mg.f:291-291] (9) 残差計算と予測 (10)アップサイクル内の近似解求解 call timer_stop(T_bench) [mg_mpi:mg.f:293-293] call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) [mg_mpi:mg.f:228-228] call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)[mg_mpi:mg.f:226-226] call mpi_reduce(t0,t,1,dp_type,[mg_mpi:mg.f:297-298] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:248-248] call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) [mg_mpi:mg.f:228-228] if( me .eq. root )then[mg_mpi:mg.f:301-365] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:249-249] if (.not.timeron) goto 999[mg_mpi:mg.f:368-368] call resid(u,v,r,n1,n2,n3,a,k) [mg_mpi :mg.f:248-248] 2024年4月 計算科学技術特論B 61 call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:256-256] do i = 1, t_last[mg_mpi:mg.f:370-372] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt)) mg_mpi :mg.f:249-249] call resid(u,v,r,n1,n2,n3,a,k)[[mg_mpi :mg.f:257-257] call MPI_Reduce(t1, tsum, t_last+2, dp_type, MPI_SUM,[mg_mpi:mg.f:376-377] call setup(n1,n2,n3,k) [mg_mpi :mg.f:258-258] call [ MPI_Reduce(t1, tming, t_last+2, dp_type, MPI_MIN,[mg_mpi call mg3P(u,v,r,a,c,n1,n2,n3,k) mg_mpi :mg.f:256-256] Kernel? Code :mg.f:378-379] < > Clear アプリケーションの構造の調査(確認) call MPI_Reduce(t1, tmaxg, t_last+2, dp_type, MPI_MAX,[mg_mpi:mg.f:380-381] call zero3(u,n1,n2,n3) [mg_mpi :mg.f:259-259] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:257-257] if (me call .eq. zran3(v,n1,n2,n3,nx(lt),ny(lt),k) 0) then[mg_mpi:mg.f:383-389] [mg_mpi:mg.f:260-260] call setup(n1,n2,n3,k)[mg_mpi :mg.f:258-258] call mpi_finalize(ierr) [mg_mpi:mg.f:395-395] call timer_stop(T_init)[mg_mpi:mg.f:262-262] subroutine show_l(z,n1,n2,n3)[show_l:mg.f:2308-2342] call zero3(u,n1,n2,n3)[mg_mpi:mg.f:259-259] do i = 1, t_last[mg_mpi:mg.f:269-271] subroutine show(z,n1,n2,n3)[show:mg.f:2386-2416] call mpi_barrier(mpi_comm_world,ierr) rprj3:制限補間 call zran3(v,n1,n2,n3,nx(lt),ny(lt),k) [mg_mpi:mg.f:260-260] [mg_mpi:mg.f:272-272] rprj3:制限補間 rprj3:制限補間 call timer_start(T_bench)[mg_mpi:mg.f:274-274] call timer_stop(T_init)[mg_mpi:mg.f:262-262] call rprj3(r(ir(k)),m1(k),m2(k),m3(k),r(ir(j)),m1(j),m2(j),m3(j),k) call rprj3(r(ir(k)),m1(k),m2(k),m3(k),r(ir(j)),m1(j),m2(j),m3(j),k) call rprj3(r(ir(k)),m1(k),m2(k),m3(k),r(ir(j)),m1(j),m2(j),m3(j),k) call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:276-276] do i = 1, t_last[mg_mpi:mg.f:269-271] call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:277-277] call mpi_barrier(mpi_comm_world,ierr) [mg_mpi do it=1,nit [mg_mpi:mg.f:272-272] :mg.f:281-288] subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k ))) subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k call :mg.f:274-274] mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:286-286] call timer_start(T_bench)[mg_mpi subroutine mg3P(u,v,r,a,c,n1,n2,n3,k)[mg3p:mg.f:570-632] (1) rprj3を調査の結果制限補間処理と判明 call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:276-276] mi(1,k),mi(2,k),mi(3,k) mi(1,k),mi(2,k),mi(3,k) mi(1,k),mi(2,k),mi(3,k) (2) ファイングリッドからコーズグリッドへの lt=lt_default=9 call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi :mg.f:277-277] lt=lt_default=9 lt=lt_default=9 do do it=1,nit[mg_mpi:mg.f:281-288] do kk==lt,1,-1 lt,1,-1 ax = 1,3 do k do = lt,1,-1 do ax = 1,3 mi(ax,k) ・・ call mg3P(u,v,r,a,c,n1,n2,n3,k) [mg_mpi:mg.f:286-286] do ax ==1,3 kk mi(ax,k) = ・・ k mip(ax,k) mi(ax,k)==・・ ・・ mip(ax,k) = ・・ subroutine mg3P(u,v,r,a,c,n1,n2,n3,k)[mg3p:mg.f:570-632] enddo mip(ax,k) = ・・ enddoenddo do k= j lt, lb+1 , -1[mg3p:mg.f:594-598] enddoenddo 補間処理を確認した (3) m1j,m2j,m3jはm1,m2,m3に結合 j j (5) その他 interp(延長補間),resid(残差計 算）,psinv（近似解求解)と確認 enddosetup call rprj3(r(ir(k)),m1(k),m2(k),m3(k),[mg.f:596-597] setup setup subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k )[rprj3:mg.f:777-869] m1,m2,m3はプロセス分散されている (6) rprj3のループ構造調査 m1,m2,m3はプロセス分散されている if (timeron) call timer_start(t_rprj3)[rprj3:mg.f:803-803] m1,m2,m3はプロセス分散されている (7) m3j,m2j,m1jの３重ループ構造 do j3=2,m3j-1[rprj3:mg.f:822-853]Br:0,St:13,FOp:24, ES2 do j2=2,m2j-1[rprj3:mg.f:825-852] do j1=2,m1j[rprj3:mg.f:829-836] do j1=2,m1j-1[rprj3:mg.f:838-850] if (timeron) call timer_stop(t_rprj3)[rprj3:mg.f:854-854] call comm3(s,m1j,m2j,m3j,j)[mg.f:858-858] 2024年4月 計算科学技術特論B 62 if( debug_vec(0) .ge. 1 )then[rprj3:mg.f:860-862] if-then-block[rprj3:mg.f:860-861] AR2:5,DAR2:5,IAR2:0,B/F:1.67 L

call resid(u(ir(k)),r(ir(k)),r(ir(k)),m1(k),m2(k),m3(k),a,k)[mg.f:618-618] do i3=1,mm3-1[interp:mg.f:994-1020]Br:0,St:11,F
call setup(n1,n2,n3,k)[mg_mpi:mg.f:224-224]
do i2=d2,mm2-1[interp:mg.f:995-1005]
subroutine resid( u,v,r,n1,n2,n3,a,k )[resid:mg.f:707-772]
Kernel?
Code
< > Clear
subroutine setup(n1,n2,n3,k)[setup:mg.f:401-565]
if [(timeron)
call timer_start(t_resid)[resid:mg.f:733-733]
call zero3(u,n1,n2,n3)
mg_mpi:mg.f:225-225]

do i1=d1,mm1-1[interp:mg.f:996-999]

アプリケーションの構造の調査(確認)

ES2
do i3=2,n3-1[resid
:mg.f:734-755]Br:0,St:7,FOp:15,
call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)
[mg_mpi
:mg.f:226-226]

do i1=1,mm1-1[interp:mg.f:1000-1004] Loop
AR2:6,DAR2:6,IAR2:0,B/F:3.20
do i2=1,mm2-1[interp:mg.f:1006-1019]

do i2=2,n2-1[resid:mg.f:735-754]
call norm2u3(v,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))
[mg_mpi:mg.f:228-228]

do i1=d1,mm1-1[interp:mg.f:1007-1011]

i1=1,n1
[resid:mg.f:736-741]
call resid(u,v,r,n1,n2,n3,a,k)do
[mg_mpi
:mg.f:248-248]

(1) psinvのループ構造調査
(1) residのループ構造調査
do i1=1,mm1-1[interp:mg.f:1012-1018]
(2)
n3,n2,n1の３重ループ構造
(2)
n3,n2,n1の３重ループ構造
if (timeron) call timer_stop(t_interp)[interp:mg.f:1023if (timeron)
call timer_stop(t_resid)
[resid:mg.f:756-756]
call mg3P(u,v,r,a,c,n1,n2,n3,k)
[mg_mpi
:mg.f:256-256]
(3) n1,n2,n3はm1,m2,m3に結合
(3) n1,n2,n3はm1,m2,m3に結合
do i1=2,n1-1[resid:mg.f:742-753]
call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))
[mg_mpi:mg.f:249-249]

call resid(u,v,r,n1,n2,n3,a,k)
[mg_mpi:mg.f:257-257]
call comm3(r,n1,n2,n3,k)
[mg.f:761-761]

call comm3_ex(u,n1,n2,n3,k)[mg.f:1025-1025]

call setup(n1,n2,n3,k)
mg_mpi:mg.f:258-258]
if([debug_vec(0)
.ge. 1 )then[resid:mg.f:763-765]

if( debug_vec(0) .ge. 1 )then[interp:mg.f:1027-1030]

call zero3(u,n1,n2,n3)
mg_mpi:mg.f:259-259]
if([debug_vec(2)
.ge. k )then[resid:mg.f:767-769]

if( debug_vec(5) .ge. k )then[interp:mg.f:1032-1035]

call zran3(v,n1,n2,n3,nx(lt),ny(lt),k)
[mg_mpi:mg.f:260-260] [mg.f:622-622]
call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)
call resid(u(ir(k)),r(ir(k)),r(ir(k)),m1(k),m2(k),m3(k),a,k)[mg.f:618call timer_stop(T_init)
[mg_mpi
:mg.f:262-262][psinv:mg.f:637-702] subroutine resid( u,v,r,n1,n2,n3,a,k )[resid:mg.f:707-772]
subroutine
psinv( r,u,n1,n2,n3,c,k)
Kernel?
do i = 1, t_last[mg_mpi
:mg.f:269-271]
if (timeron) call timer_start(t_psinv)[psinv:mg.f:663-663]

if (timeron) call timer_start(t_resid)[resid:mg.f:733-73

call mpi_barrier(mpi_comm_world,ierr)[mg_mpi:mg.f:272-272]
Loo
do i3=2,n3-1[psinv:mg.f:664-685]Br:0,St:7,FOp:16, ES2
AR2:6,DAR2:6,IAR2:0,B/F:3.00
do i3=2,n3-1
[resid:mg.f:734-755]Br:0,St:7,FOp:15, ES2
call timer_start(T_bench)[mg_mpi:mg.f:274-274]
do i2=2,n2-1[psinv:mg.f:665-684]
do i2=2,n2-1[resid:mg.f:735-754]
call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:276-276]
do i1=1,n1[psinv:mg.f:666-671]
do i1=1,n1[resid:mg.f:736-741]
call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:277-277]
do i1=2,n1-1[psinv:mg.f:672-683]
do i1=2,n1-1[resid:mg.f:742-753]
do it=1,nit[mg_mpi:mg.f:281-288]
if (timeron) call timer_stop(t_psinv)[psinv:mg.f:686-686]
if (timeron) call timer_stop(t_resid)[resid:mg.f:756-756
call mg3P(u,v,r,a,c,n1,n2,n3,k)[mg_mpi:mg.f:286-286]
call comm3(u,n1,n2,n3,k)[mg.f:691-691]
call comm3(r,n1,n2,n3,k)[mg.f:761-761]
subroutine mg3P(u,v,r,a,c,n1,n2,n3,k)[mg3p:mg.f:570-632]
if( debug_vec(0) .ge. 1 )then[psinv:mg.f:693-695]
if( debug_vec(0) .ge. 1 )then[resid:mg.f:763-765]
do k= lt, lb+1 , -1[mg3p:mg.f:594-598]
if-then-block[psinv:mg.f:693-694]
if( debug_vec(2) .ge. k )then[resid:mg.f:767-769]
call rprj3(r(ir(k)),m1(k),m2(k),m3(k),[mg.f:596-597]
call rep_nrm(u,n1,n2,n3,' psinv',k)[mg.f:694-694]
call :mg.f:777-869]
psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)[mg.f:622-622]
subroutine rprj3( r,m1k,m2k,m3k,s,m1j,m2j,m3j,k )[rprj3
if( debug_vec(3) .ge. k )then[psinv:mg.f:697-699]
subroutine psinv( r,u,n1,n2,n3,c,k)[psinv:mg.f:637-702]
if (timeron) call timer_start(t_rprj3)[rprj3:mg.f:803-803]
call interp(u(ir(j)),m1(j),m2(j),m3(j),u,n1,n2,n3,k)
[mg.f:627-627]
Loo
do j3=2,m3j-1[rprj3:mg.f:822-853]Br:0,St:13,FOp:24,if ES2
AR2:5,DAR2:5,IAR2:0,B/F:1.67
(timeron)
call timer_start(t_psinv)[psinv:mg.f:663-66

call resid(u,v,r,n1,n2,n3,a,k)
[mg.f:628-628]
do j2=2,m2j-1
[rprj3:mg.f:825-852]
2024年4月 計算科学技術特論B

call psinv(r,u,n1,n2,n3,c,k)
[mg.f:629-629]
do j1=2,m1j
[rprj3:mg.f:829-836]

do i3=2,n3-1[psinv:mg.f:664-685]Br:0,St:7,FOp:16, ES2

63

do i2=2,n2-1[psinv:mg.f:665-684]

call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi
:mg.f:287-287]
do j1=2,m1j-1
[rprj3:mg.f:838-850]

do i1=1,n1[psinv:mg.f:666-671]

アプリケーションの構造の調査(確認)

call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))
[mg_mpi:mg.f:291-291]
if (timeron) call timer_stop(t_rprj3)
[rprj3:mg.f:854-854]

call timer_stop(T_bench)call
[mg_mpi
:mg.f:293-293][mg.f:858-858]
comm3(s,m1j,m2j,m3j,j)

do i1=2,n1-1[psinv:mg.f:672-683]

if (timeron) call timer_stop(t_psinv)[psinv:mg.f:686-68

if( debug_vec(0)
.ge. 1 )then[rprj3:mg.f:860-862]
call mpi_reduce(t0,t,1,dp_type,
[mg_mpi:mg.f:297-298]

call comm3(u,n1,n2,n3,k)[mg.f:691-691]

(1) interpのループ構造調査
if( debug_vec(0) .ge. 1 )then[psinv:mg.f:693-695]
(2) mm3,mm2,mm1の３重ループ構造
call rep_nrm(s,m1j,m2j,m3j,' rprj3',k-1)[mg.f:861-861]
if (.not.timeron) goto 999[mg_mpi
:mg.f:368-368]
if-then-block[psinv:mg.f:693-694]
mm1,mm2,mm3はm1,m2,m3に結合
if( (3)
debug_vec(4)
.ge. k )then[rprj3:mg.f:864-866]
do i = 1, t_last[mg_mpi:mg.f:370-372]

if-then-block[rprj3:mg.f:860-861]
if( me .eq. root )then[mg_mpi:mg.f:301-365]

call zero3(u(ir(k)),m1(k),m2(k),m3(k))
[mg.f:604-604]
call MPI_Reduce(t1,
tsum, t_last+2, dp_type, MPI_SUM,
[mg_mpi:mg.f:376-377]

call rep_nrm(u,n1,n2,n3,' psinv',k)[mg.f:694-694]

if( debug_vec(3) .ge. k )then[psinv:mg.f:697-699]
call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)
[mg.f:605-605]
call MPI_Reduce(t1,
tming, t_last+2, dp_type, MPI_MIN,[mg_mpi
:mg.f:378-379]
call interp(u(ir(j)),m1(j),m2(j),m3(j),u,n1,n2,n3,k)[mg.f:627-627]
do k = lb+1, lt-1[mg3p:mg.f:607-623]
call MPI_Reduce(t1, tmaxg, t_last+2, dp_type, MPI_MAX,[mg_mpi:mg.f:380-381]
call resid(u,v,r,n1,n2,n3,a,k)[mg.f:628-628]
call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:612-612]
if (me .eq. 0) then[mg_mpi:mg.f:383-389]
call interp(u(ir(j)),m1(j),m2(j),m3(j),[mg.f:613-614]call psinv(r,u,n1,n2,n3,c,k)[mg.f:629-629]
call mpi_finalize(ierr)[mg_mpi:mg.f:395-395]
call resid(u,v,r,n1,n2,n3,a,k)
[mg_mpi:mg.f:287-287]
subroutine interp( z,mm1,mm2,mm3,u,n1,n2,n3,k
)[interp:mg.f:874-1038]
outine show_l(z,n1,n2,n3)[show_l:mg.f:2308-2342]
if (timeron) call timer_start(t_interp)
[interp:mg.f:905-905]
call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))
[mg_mpi:mg.f:291-291]
if( n1 .ne. 3 .and. n2 .ne. 3 .and. call
n3 .ne.
3 ) then[interp:mg.f:906-1022]
timer_stop(T_bench)
[mg_mpi:mg.f:293-293]
if-then-block[interp:mg.f:906-942]
call mpi_reduce(t0,t,1,dp_type,[mg_mpi:mg.f:297-298]
ES2
do i3=1,mm3-1[interp:mg.f:908-942]Br:0,St:18,FOp:23,
AR2:8,DAR2:8,IAR2:0,B/F:2.78
if( me .eq. root )then[mg_mpi:mg.f:301-365]

do i2=1,mm2-1[interp:mg.f:909-941]

if (.not.timeron) goto 999[mg_mpi:mg.f:368-368]

do i1=1,mm1[interp:mg.f:911-915]

do i = 1, t_last[mg_mpi:mg.f:370-372]

do i1=1,mm1-1[interp:mg.f:917-922]

(1) 主要４サブルーチンとも処理量は
m1,m2,m3の大きさに依存している．

call MPI_Reduce(t1, tsum, t_last+2, dp_type, MPI_SUM,[mg_mpi:mg.f:376-377]

do i1=1,mm1-1[interp:mg.f:923-928]

call:mg.f:929-934]
MPI_Reduce(t1, tming, t_last+2, dp_type, MPI_MIN,[mg_mpi:mg.f:378-379]
do i1=1,mm1-1[interp

(2) m1,m2,m3の分割はどうなっているか？

call:mg.f:935-940]
MPI_Reduce(t1, tmaxg, t_last+2, dp_type, MPI_MAX,[mg_mpi:mg.f:380-381]
do i1=1,mm1-1[interp
else-block[interp:mg.f:944-1020]
if (me .eq. 0) then[mg_mpi:mg.f:383-389]
ES2
do i3=d3,mm3-1[interp:mg.f:970-992]Br:0,St:11,FOp:12,
AR2:7,DAR2:7,IAR2:0,B/F:4.67
call mpi_finalize(ierr)[mg_mpi:mg.f:395-395]
2024年4月

計算科学技術特論B

64
do i2=d2,mm2-1[interp:mg.f:971-980]

subroutine show_l(z,n1,n2,n3)[show_l:mg.f:2308-2342]

do i1=d1,mm1-1[interp:mg.f:972-975]

データ分割の調査 (1) setup内の配列：miを調べた結 果，プロセス分割されている事 が判明 2024年4月 計算科学技術特論B 65 スケーラビリティの調査 (1) ４プロセスの測定．ウィークスケーリング測定 主要処理のelaps時間 全体で17sec Performance monitor event ************************************************************************************* Application - performance monitors ************************************************************************************* Elapsed(s) MFLOPS MFLOPS/PEAK(%) MIPS MIPS/PEAK(%) -----------------------------------------------------------------------16.6556 10361.0515 2.0236 18605.9895 7.2680 -----------------------------------------------------------------------16.6556 2590.2629 2.0236 4651.6009 7.2681 16.6052 2598.1206 2.0298 4668.2023 7.2941 16.5407 2608.2634 2.0377 4683.4974 7.3180 16.5213 2611.3249 2.0401 4686.9329 7.3233 Mem throughput Mem throughput Elapsed(s) _chip(GB/S) /PEAK(%) SIMD(%) --------------------------------------------------------16.6556 32.3398 12.6327 48.7947 --------------------------------------------------------16.6556 8.0885 12.6382 48.7936 16.6052 8.1070 12.6672 48.7676 16.5407 8.1357 12.7120 48.7981 16.5213 8.1551 12.7423 48.8196 2024年4月 計算科学技術特論B Application Process Process Process Process Application Process Process Process Process 66 0 2 1 3 0 2 1 3 プロセスのロード インバランスは 発生していない ピーク性能比２%

スケーラビリティの調査 (1) 32プロセスの測定．ウィークスケーリング測定 主要処理のelaps時間 全体で40sec Performance monitor event ************************************************************************************* Application - performance monitors ************************************************************************************* Elapsed(s) MFLOPS MFLOPS/PEAK(%) MIPS MIPS/PEAK(%) -----------------------------------------------------------------------39.3654 83044.5238 2.0275 149003.2766 7.2756 -----------------------------------------------------------------------39.3654 2595.1412 2.0275 4654.6126 7.2728 39.3575 2595.6627 2.0279 4650.9112 7.2670 39.3546 2595.8520 2.0280 4656.6897 7.2761 : Mem throughput Mem throughput Elapsed(s) _chip(GB/S) /PEAK(%) SIMD(%) --------------------------------------------------------39.3654 258.7562 12.6346 49.3192 --------------------------------------------------------39.3654 8.0773 12.6208 49.3303 39.3575 8.0573 12.5895 49.3647 39.3546 8.0544 12.5850 49.3071 : 2024年4月 計算科学技術特論B プロセスのロード インバランスは 発生していない ピーク性能比２% Application Process Process Process 0 7 8 Application Process Process Process 0 7 8 67 スケーラビリティの調査 ウィークスケーリングで４プロセスから32プロセスで全体で17sec から40secへ実行時間が増大．約2.5倍の増大． 主要４サブルーチンも同様の傾向． しかしピーク性能比は、２％で変わらず． １プロセスの演算量が同じなら実行時間が延びる事でピーク性能 比も落ちるはず？ 調査の結果ウィークスケール測定であったが32プロセスでは全 体のイタレーション回数が2.5倍に設定されていた事が判明． スケーラビリティとしては良好である事が分かった． またロードインバランスも良好であることが分かった． 並列化の観点では問題なし． 2024年4月 計算科学技術特論B 68

高並列特性評価(結果) 実行時間 物理的 ・スケーラ 処理内容 ビリティ do k= lt, lb+1 , -1[mg3p :mg.f:594-598] 演算・通 信特性 演算・通 信見積り カー ネル call rprj3(r(ir(k)),m1(k),m2(k),m3(k),[mg.f:596-597] 良好 制限補間 完全並列 Nに比例 call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:604-604] ○ call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k)[mg.f:605-605] do k = lb+1, lt-1[mg3p:mg.f:607-623] call zero3(u(ir(k)),m1(k),m2(k),m3(k))[mg.f:612-612] 完全並列 Nに比例 call interp(u(ir(j)),m1(j),m2(j),m3(j), 良好 延長補間[mg.f:613-614] ○ call resid(u(ir(k)),r(ir(k)),r(ir(k)),m1(k),m2(k),m3(k),a,k) [mg.f:618-618] 残差計算 完全並列 Nに比例 良好 ○ call psinv(r(ir(k)),u(ir(k)),m1(k),m2(k),m3(k),c,k) 簡易求解 完全並列 [mg.f:622-622] Nに比例 良好 ○ call interp(u(ir(j)),m1(j),m2(j),m3(j),u,n1,n2,n3,k)[mg.f:627-627] call resid(u,v,r,n1,n2,n3,a,k)[mg.f:628-628] call psinv(r,u,n1,n2,n3,c,k)[mg.f:629-629] call resid(u,v,r,n1,n2,n3,a,k)[mg_mpi:mg.f:287-287] 2024年4月 計算科学技術特論B 69 call norm2u3(r,n1,n2,n3,rnm2,rnmu,nx(lt),ny(lt),nz(lt))[mg_mpi:mg.f:291-291] まとめ call timer_stop(T_bench)[mg_mpi:mg.f:293-293] call mpi_reduce(t0,t,1,dp_type,[mg_mpi:mg.f:297-298] if( me .eq. root )then[mg_mpi:mg.f:301-365] if (.not.timeron) goto 999[mg_mpi:mg.f:368-368] • 反復法 do i = 1, t_last[mg_mpi:mg.f:370-372] • 並列化できるアルゴリズムとして反復法について説明しました． call MPI_Reduce(t1, tsum, t_last+2, dp_type, MPI_SUM,[mg_mpi:mg.f:376-377] call MPI_Reduce(t1, tming, t_last+2, dp_type, MPI_MIN,[mg_mpi:mg.f:378-379] • アプリケーションの性能最適化 call MPI_Reduce(t1, tmaxg, t_last+2, dp_type, MPI_MAX,[mg_mpi:mg.f:380-381] 性能最適化の流れと高並列化のための分析手法を説明しました． if (me• .eq. 0) then[mg_mpi:mg.f:383-389] call mpi_finalize(ierr)[mg_mpi:mg.f:395-395] • 並列性能のボトルネック subroutine show_l(z,n1,n2,n3)[show_l:mg.f:2308-2342] subroutine show(z,n1,n2,n3)[show:mg.f:2386-2416] • 並列性能のボトルネック要因について説明し例を示しました • 簡単な実例 • NPB MGを例に高並列化のための分析手法の実例を示しました． 2024年4月 計算科学技術特論B 70

= A - - F E D −4 1 0 0 0 1 −4 1 0 0 A= 0 1 −4 1 0 0 0 1 −4 1 0 0 0 1 −4 = D-E 2024年4月 2022年4月 計算科学技術特論B 71 2024年4月 計算科学技術特論B 72