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

計算科学特論B 第5回 アプリケーションの性能最適化の実例２ ジャパンメディカルデバイス株式会社 熊畑 清 2024.5.16(⽊)

講義の概要 • 4⽉11⽇ スーパーコンピュータとアプリケーションの性能 • 4⽉18⽇ アプリケーションの性能最適化１(⾼並列性能最適化) • 4⽉25⽇ アプリケーションの性能最適化の実例1とCPU単体性能 • 5⽉9⽇ アプリケーションの性能最適化2（CPU単体性能最適化） • 5⽉16⽇ アプリケーションの性能最適化の実例2 2

本⽇の内容 • • 流体解析ソルバーでの例 • • キャッシュ利⽤効率の改善・スレッド並列化の⽅法の⼯夫 ストアのシーケンシャル化 • ストリーム削減・セクタキャッシュ・プリフェッチ 乱流燃焼コードでの例 • • インライン展開のやり過ぎに注意 ⼩回転ループはさけたい • • OpenMPのスケジュール調整 重複した処理の回避 • 配列初期化時の性能劣化 • その他Tips • なるべく⻑いループでSIMDに • Z-Fillの効果 • ifを含むループはリストベクトル化 • 多次元配列の配列初期化 3

流体解析ソルバー キャッシュ利⽤効率の改善 スレッド並列化の⽅法の⼯夫

チューニング対象カーネル 要素圧力モードでの処理時間内訳，上位10件 速度 1.89 2.10 2.95 4.83 14.11 全体 #$%#$ 単位( 秒 8.32 13.04 callap._PRL_2_ calaxc_ fild3x._PRL_1_ callap_ vel3d1._PRL_16_ int3dx_ callap._PRL_5_ vel3d1._PRL_8_ vel3d1_ nodlex_ 速度・圧⼒ 圧⼒ 速度 速度 速度・圧⼒ 要素圧⼒モード 速度・圧⼒ 節点圧⼒モード • 10万個の4面体要素＋α • 14.11s(30%)が4面体要素の勾配計算 • 13.04s(27%)が疎行列ベクトル積 • 8.32s（17%)が4面体要素の発散計算 • Callap._PRL_2_のチューニングを実施 5

チューニング対象カーネル DO ICOLOR=1,NCOLOR(1) IES=LLOOP(ICOLOR ,1)+1 IEE=LLOOP(ICOLOR+1,1) DO IE=IES,IEE・・・・・・・ここで並列 IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO 4角形要素について∇pの有限要素近似 を計算するループ ¶N j ¶p Ñp = @å pe ¶xi j =1 ¶xi 型 配列名とサイズ INTEGER*4 NODE(9,NE) REAL*4 S(NE) DNX(9,NE), DNY(9,NE), REAL*4 DNZ(9,NE) REAL*4 FX(NP), FY(NP), FZ(NP) 内容 節点リスト 圧力 形状関数の導関数 圧力勾配ベクトル 6

チューニング対象カーネル DO ICOLOR=1,NCOLOR(1) IES=LLOOP(ICOLOR ,1)+1 IEE=LLOOP(ICOLOR+1,1) DO IE=IES,IEE・・・・・・・ここで並列 IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO 4角形要素について∇pの有限要素近似 を計算するループ ¶N j ¶p Ñp = @å pe ¶xi j =1 ¶xi FX(1) FX(2) FX(3) FX(4) S(1)*DNX(1,1) S(1)*DNX(2,1) S(2)*DNX(1,2) S(2)*DNX(2,2) S(3)*DNX(1,3) S(3)*DNX(2,3) S(1)*DNX(4,1) S(1)*DNX(3,1) S(2)*DNX(4,2) S(2)*DNX(3,2) S(3)*DNX(4,3) S(3)*DNX(3,3) FX(5) FX(6) FX(7) FX(8) 7

チューニング対象カーネル DO ICOLOR=1,NCOLOR(1) IES=LLOOP(ICOLOR ,1)+1 IEE=LLOOP(ICOLOR+1,1) DO IE=IES,IEE・・・・・・・ここで並列 IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO カラーリングによりデータ依存が⽣ずる 要素を別のグループ(カラー)に分ける 1 5 2 6 1 5 9 13 10 14 9 13 10 14 3 7 3 7 4 8 11 15 12 16 2 4 6 8 11 15 12 16 要素(メッ シュ) 配列 LLOOP/NODE/DNX〜Zメモリ上配置 1 2 3 4 カラー1 5 6 7 8 ・・・ カラー2 ・・・ 9 10 11 12 13 14 15 16 カラー3 カラー4 8

性能評価 要求バイト数の算出 配列FX,FY,FZ • 配列NODEの値に依存するリストアクセス故にシーケンシャルではない • 再利用性がある • リストによって変化してしまい一概には算出できない DO IE=IES, IEE IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) FX(1) FX(2) FX(4) FX(3) リストアクセス FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) ENDDO S(4)*DNX(1,4) S(4)*DNX(2,4) FX(5) S(N)*DNX(2,N) S(2)*DNX(1,2) S(4)*DNX(3,4) S(N)*DNX(4,N) FX(6) S(2)*DNX(2,2) IE=2 IE=N IE=1 S(4)*DNX(4,4) S(N)*DNX(1,N) S(N)*DNX(3,N) S(2)*DNX(4,2) FX(7) S(2)*DNX(3,2) FX(8) 9

性能評価 要求バイト数の算出 配列FX,FY,FZ • 配列NODEの値に依存するリストアクセス故にシーケンシャルではない • 再利用性がある • リストによって変化してしまい一概には算出できない IE=IE+1 YES DO IE=IES, IEE 終了か IP1=NODE(1,IE) IP2=NODE(2,IE) NO IP3=NODE(3,IE) IP1=NODE(1,IE) IP4=NODE(4,IE) FX(IP1) 最後 リストアクセス IE=1 IE=2 IE=3 FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) キャッシュ 4 1,3 10 ENDDO 再利⽤性 FX(100) IE=4 IE=5 ・・・・ IE=ME 100 3 2 Elem 1 2 10 25 n 5 1 IEごとに触る 4 Elem 配列FXの成分 2 3 51 全てこの範囲で収まる のが理想的 従って理論性能の算出に当たり，メモリプレッシャーはないと仮定する 10

性能評価 要求バイト数の算出 配列S • 再内のループ変数IEに従い連続アクセス • 再利用性はない • 1回目のアクセスでキャッシュミスし1ライン分128Byteがメモリ から転送される DO IE=IES,IEE ・・・ SWRK=S(IE) ・・・ ENDDO • ループ1回転で4Byte、128Byteでは32回転を賄える • メモリプレッシャーは32回転につき128Byte ・・・ 17 18 19 6 20 7 21 8 22 9 23 10 24 11 25 12 26 13 27 14 28 15 29 16 30 ・・・ 31 32 ↑ S(1)でキャッシュミス ↑ 1 2 3 4 5 S(32)まで1ラインで賄える 1ライン分128Byte＝32個 11

性能評価 要求バイト数の算出 配列NODE INTEGER*4 NODE(9,NE) • 再内のループ変数IEに従い規則的にアクセス • 再利用性はない • 1回目のアクセスでキャッシュミスし1ライン分128Byteがメモリ から転送される • 第1次元のサイズが9でありループ1回転で36Byte、128Byte では3.56回転を賄える • 32回転で9回ミスしメモリプレッシャーは32回転につき 1152Byte NODE(1,1) ↑ 1 2 3 ・・・ 17 4 18 5 19 ↓ 6 20 7 21 NODE(1,3) 8 22 9 23 NODE(1,2) ↑ 10 11 12 24 25 26 13 27 14 28 ↓ DO IE=IES,IEE IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) ENDDO 15 29 16 30 NODE(1,4) ・・・ 31 32 12

性能評価 要求バイト数の算出 配列DNX, DNY, DNZ INTEGER*4 NODE(9,NE) • 第1次元のサイズが9であり配列NODEと同じ性質 • メモリプレッシャーは32回転につきそれぞれ1152Byte DO IE=IES,IEE ・・DNX(1,IE) ・・DNX(2,IE) ・・DNX(3,IE) ・・DNX(4,IE) ENDDO 使われないがメモリ帯域を使う分を考慮する必要 DNX(1,1) ↑ 1 2 3 4 5 6 7 8 DNX(1,2) ↑ 10 11 12 9 DNX(5～9,1)は使われない ・・・ 17 18 19 ↓ 20 21 NODE(1,3) 22 23 24 13 14 15 16 ・・・ DNX(5～9,2)は使われない 25 26 27 DNX(5～9,3)は使われない 28 ↓ 29 30 NODE(1,4) 31 32 13

性能評価 • 要求バイト 128＋1152＋1152×3＝4736Byte 5.0 • 浮動小数点演算数 32回転×24FLOP＝768Flop 3.0 命令実行 2.0 待ち 1.0 メモリアクセス関連待 ち • B/F値は 4736/768=6.17 • STREAMベンチマークによる「京」の 実効メモリバンド幅は46.6GB/s • 「京」の実効B/F値は46.6/128=0.36 • メモリバンド幅がネックとなり実効上 の性能上限は128GFLOPSの 0.36/6.17=5.83% 計算時間内訳 Time [s] 4.0 0.0 1 2 3 4 5 6 7 8 スレッド • 82万個の4面体要素 • マルチスレッドで1.6% • メモリスループットは10.29GB/s • メモリアクセス関連待ちが70%程度 • L1Dキャッシュミス率21% • キャッシュ利用効率が悪い 14

キャッシュの利⽤効率改善 • 配列FX, FY, FZへのアクセスには再利用性がある • 常にオンキャッシュであると仮定(理想的) • 実際の入力データでは再利用性はあるものの，幾何的に近い節点に，近い番号がつけ られているとは限らない 節点 1 18 1 4 要素１ 2 75 3 40 参照 #$$%&'' FX( 1) + FX(75) ,-. FX(40) /0. FX(18) .1 ()*& + , / + 1=NODE(1,IE) 18=NODE(2,IE) 40=NODE(3,IE) 互いに遠い節点を参照(跳び) 75=NODE(4,IE) キャッシュの利⽤効効率が悪い ・・・ FX( 1)=FX( 1)-SWRK*DNX(1,IE) FX(18)=FX(18)-SWRK*DNX(2,IE) FX(40)=FX(40)-SWRK*DNX(3,IE) FX(75)=FX(75)-SWRK*DNX(4,IE) 幾何的に近い節点には近い番号をつける必要 15

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ① ③ Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ⑨ ⑥ Z ⑪ Y 椃 Axial-Aligned Bounding Box （AABB) 16

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ① ① Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ⑨ ⑥ Z ⑪ Y 椃 Axial-Aligned Bounding Box （AABB) 17

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ① ① Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ⑨ ② Z ⑪ Y 椃 Axial-Aligned Bounding Box （AABB) 18

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ① ① Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ③ ② Z ⑪ Y 椃 Axial-Aligned Bounding Box （AABB) 19

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ④ ① Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ③ ② Z ⑪ Y 椃 Axial-Aligned Bounding Box （AABB) 20

キャッシュの利⽤効率改善 • 幾何的に近い節点には近い番号をつける • 行列・ベクトル席の前処理として，バンド幅の縮小や，直接法を用いる際のFill-inの抑 制を目的として行や列の並びを変更する手法 • Minimum Degree法, Nested Dissection法, グラフ理論を応用した方法など • 行列を作らずに要素毎に計算する実装なため行や列の並べ替え操作は，節点の番 号を入れ替えることに相当 ( 外 ④ ① Object 10 11 5 棙 4 棖 1 棑 2 棒 3 棓 ' 中 ③ ② Z ⑤ Y 椃 Axial-Aligned Bounding Box （AABB) 21

キャッシュの利⽤効率改善 • 幾何的に近い節点に近い番号がつけられた • 幾何的に近い節点の情報が，メモリ上でも近い位置になるようにメモリ上でソートして おく 節点○の配置 ② 直⽅体1 ① ③ ⑤ ⑧ ④ ⑨ ⑥ ⑦ ⑪ ⑩ リオーダリング ④ ① ⑤ ② ③ ⑥ ⑪ ⑩ ⑦ ⑧ ⑨ 直⽅体2 index 節点の情報を保持する配列の並び 1 2 3 4 5 6 7 8 9 10 11 節点情報 ① ② ③ ④ ⑤ ⑥ ⑦ ⑧ ⑨ ⑩ ⑪ index 1 2 3 4 5 6 7 8 9 10 11 節点情報 ④ ⑩ ① ⑦ ⑤ ② ⑧ ⑨ ③ ⑪ ⑥ index 1 2 3 4 5 6 7 8 9 10 11 節点情報 ① ② ③ ④ ⑤ ⑥ ⑦ ⑧ ⑨ ⑩ ⑪ 幾何的に近い節点で もメモリ上は遠い 番号の付け直しと メモリ上でのソート 幾何的に近い節点は メモリ上でも近い 22

Appendix) リオーダーの効果 同じアプリの別部分 DO I=1,N IDX=LIST(I) V=ARRAY(IDX) ・・・ ENDDO 図は、横軸がループの回転、縦軸が、そのループで アクセスする配列インデックスを⽰したもの e r o f Be ４面体要素 ６面体要素 ※一部抜粋 23

Appendix) リオーダーの効果 同じアプリの別部分 DO I=1,N IDX=LIST(I) V=ARRAY(IDX) ・・・ ENDDO 図は、横軸がループの回転、縦軸が、そのループで アクセスする配列インデックスを⽰したもの A r e ft ４面体要素 ６面体要素 ICRS=1～5M 24

Appendix) リオーダーの効果 同じアプリの別部分 After(6⾯対) スレッド4 スレッド3 スレッド2 スレッド1 スレッド0 Before(6⾯体) 25

スレッド並列化の⽅法の⼯夫 • 要素のカラーリング領域全体で行っているので，最内ループで中には幾何的に遠い 要素が出現する • 参照する節点がメモリ上でも遠い位置あり不利 メッシュ カラーリング 要素1個 26

スレッド並列化の⽅法の⼯夫 • 要素のカラーリング領域全体で行っているので，最内ループで中には幾何的に遠い 要素が出現する • 参照する節点がメモリ上でも遠い位置あり不利 最内ループ1単位分 カラーリング 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 27

スレッド並列化の⽅法の⼯夫 • 要素のカラーリング領域全体で行っているので，最内ループで中には幾何的に遠い 要素が出現する • 参照する節点がメモリ上でも遠い位置あり不利 最内ループ1単位分 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 メモリ 上 ・ ・ ・ キャッシュから外れる可能性 ・ ・ ・ カラー内に含まれる要素を，幾何的に近いものにする必要 28

スレッド並列化の⽅法の⼯夫 • カラー内に含まれる要素を，幾何的に 近いものにする • ブロックで動くループが追加 1ノード分の 計算領域 ブロック化 □ はブロック1個 R B R G Y G R B R カラーリング □ は要素1個 DO IBLOCK=1,NBLOCK(1) DO ICOLOR=1,NCOLOR(IBLOCK,1) IES=LLOOP(IBLOCK,ICOLOR ,1)+1 IEE=LLOOP(IBLOCK,ICOLOR+1,1) DO IE=IES,IEE・・・・・・・ここで並列 IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO ENDDO 29

スレッド並列化の⽅法の⼯夫 ブロックの大きさ(含む要素数)について検討 1要素当たりの必要メモリ量 ケースNo IP1=NODE(1,IE) NODE(1～9, IE) IP2=NODE(2,IE) 36バイト IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) S(IE)→4バイト FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) DNX(1-9,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) 36バイト FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) DNZ(1-9,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) 36バイト FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) FX(IP1～IP4) FY(IP1～IP4) FZ(IP1～IP4) 48バイト 196バイト/要素 on L1 on L2 1 640 1291 247.1 30.9 ○ ○ 2 320 2583 494.4 61.8 × ○ 3 160 5166 988.8 123.6 × ○ 4 80 10333 1977.8 247.2 × ○ 5 40 20666 3955.6 494.5 × ○ 6 20 41332 7911.2 988.9 × × 7 10 82665 15822.6 1977.8 × × 8 5 165331 31645.4 3955.7 × × 9 3 275552 52742.4 6592.8 × × 10 2 413328 79113.6 9889.2 × × 11 1 826656 158227.1 19778.4 × × 1000000 100000 アクセスデータ量 [KB] FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) DNY(1-9,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) 36バイト FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) ブロック 平均要素数/ アクセスデー 1スレッド当り 分割数 ブロック タ量 [KB] メモリ [KB] 10000 1000 100 ブロック当たりに アクセスするデータ量 10 1 1000 10000 100000 1000000 平均要素数/ブロック (対数表示) 30

スレッド並列化の⽅法の⼯夫 ブロックの大きさ(含む要素数)について検討 2.00% 40000 ピーク性能比 % インデックスの飛び ピーク性能比 30000 20000 10000 0 1000 10000 100000 平均要素数/ブロック (対数表示) 1 1000000 1.50% 1.00% 0.50% 0.00% 1000 10000 100000 1000000 平均要素数/ブロック (対数表示) 301 アクセスの跳び =301-1=300 100 88 31

スレッド並列化の⽅法の⼯夫 ブロックの大きさ(含む要素数)について検討 ブロック 平均要素数/ 平均カラー数 平均要素数/ 平均最内回 分割数 ブロック /ブロック カラー 転数/スレッド 1 640 1291 37.4 34.0 4 2 320 2583 38.7 66.0 8 3 160 5166 39.7 130.0 16 4 80 10333 41.0 252.0 32 5 40 20666 41.5 498.0 62 6 20 41332 42.1 981.0 123 7 10 82665 42.2 1958.0 245 8 5 165331 42.6 3881.0 485 9 3 275552 43.7 6310.0 789 10 2 413328 43.5 9501.0 1188 11 1 826656 44.0 18787.0 2348 2500 平均最内回転数/スレッド ケースNo 2000 最内ループ平均回転数/スレッド 1500 1000 最内ループ回転数が⼩さ くなりすぎた 回転数不⾜の回避が必要 500 0 1000 10000 100000 1000000 平均要素数/ブロック (対数表示) 32

スレッド並列化の⽅法の⼯夫 • カラー内に含まれる要素を，幾何的に 近いものにする • ブロックとカラーの関係を反転 1ノード分の 計算領域 ブロック化 □ はブロック1個 カラーリング DO ICOLOR=1,NCOLOR(1) DO IBLOCK=1,NBLOCK(IBLOCK,1) IES=LLOOP(IBLOCK,ICOLOR ,1)+1 IEE=LLOOP(IBLOCK,ICOLOR+1,1) DO IE=IES,IEE IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO ENDDO 33

スレッド並列化の⽅法の⼯夫 DO ICOLOR=1,NCOLOR(1) DO IBLOCK=1,NBLOCK(IBLOCK,1) ここで並列 IES=LLOOP(IBLOCK,ICOLOR ,1)+1 IEE=LLOOP(IBLOCK,ICOLOR+1,1) DO IE=IES,IEE IP1=NODE(1,IE) IP2=NODE(2,IE) IP3=NODE(3,IE) IP4=NODE(4,IE) SWRK=S(IE) FX(IP1)=FX(IP1)-SWRK*DNX(1,IE) FX(IP2)=FX(IP2)-SWRK*DNX(2,IE) FX(IP3)=FX(IP3)-SWRK*DNX(3,IE) FX(IP4)=FX(IP4)-SWRK*DNX(4,IE) • カラー内に含まれる要素を，幾何的に 近いものにする • ブロックとカラーの関係を反転 1ノード分の 計算領域 ブロック化 □ はブロック1個 ブロックの 中には 複数の要素 ブロックを カラーリング FY(IP1)=FY(IP1)-SWRK*DNY(1,IE) FY(IP2)=FY(IP2)-SWRK*DNY(2,IE) FY(IP3)=FY(IP3)-SWRK*DNY(3,IE) FY(IP4)=FY(IP4)-SWRK*DNY(4,IE) FZ(IP1)=FZ(IP1)-SWRK*DNZ(1,IE) FZ(IP2)=FZ(IP2)-SWRK*DNZ(2,IE) FZ(IP3)=FZ(IP3)-SWRK*DNZ(3,IE) FZ(IP4)=FZ(IP4)-SWRK*DNZ(4,IE) ENDDO ENDDO ENDDO 34

スレッド並列化の⽅法の⼯夫 4.00% 最内ループ 平均回転数 100000 10000 1000 Before 100 After 10 ピーク性能⽐ 3.50% ピーク性能比 % 平均最内回転数/スレッド (対数表示) 1000000 3.00% 2.50% 2.00% 1.50% 1.00% 0.50% 1 1000 10000 100000 1000000 平均要素数/ブロック (対数表示) • 最内ループの回転数が38～44倍程度ま で増大 • ブロックサイズ小の側で性能向上 • 最良の結果は平均最内回転数/スレッド が323回転 0.00% 1000 10000 100000 1000000 平均要素数/ブロック (対数表示) • 平均最内回転数/スレッド=323回転 • ピーク性能比3.78% • L1Dミス率4％ • メモリスループット32GB/s 35

スレッド並列化の⽅法の⼯夫 ここまでのチューニングにより オリジナル 改善後 理想値 ピーク性能比(%) 1.6 3.8 5.83 メモリスループット(GB/s) 10.3 32 46 L1Dキャッシュミス率(%) 21.0 3.78 3.13 もう少し改善の余地があるのでは？ 4バイト変数の配列S(N) 2 ・・・ 3 17 4 18 5 19 6 20 7 21 8 22 9 23 10 24 11 25 12 26 13 27 14 28 15 29 16 30 ・・・ 31 32 ↑ 1 S(32)まで1ラインで賄える 1ライン分128Byte＝32個 32回に1回ミスするので1/32=3.125%が理論的なミス率 (ただし全てシーケンシャルな場合なので注意) 36

スレッド並列化の⽅法の⼯夫 配列融合によるストリーム削減 オリジナル FX(IP1)=FX(IP1)+SWRK*DNX(1,IE) FY(IP1)=FY(IP1)+SWRK*DNY(1,IE) FZ(IP1)=FZ(IP1)+SWRK*DNZ(1,IE) パターン4 FXYZ(1,IP1)=FXYZ(1,IP1)+SWRK*DNXYZ(1,1,IE) FXYZ(2,IP1)=FXYZ(2,IP1)+SWRK*DNXYZ(2,1,IE) FXYZ(3,IP1)=FXYZ(3,IP1)+SWRK*DNXYZ(3,1,IE) パターン1 パターン2 パターン3 パターン4 オリジナル 内容 パターンNo FXYZ(1,IP1)=FXYZ(1,IP1)-SWRK*DNXYZ(1,1,IE) FX(1,IP1)=FX(1,IP1)-SWRK*DNXYZ(1,IE) FXYZ(1,IP1)=FXYZ(1,IP1)-SWRK*DNX(1,IE) FX(1,IP1)=FX(1,IP1)-SWRK*DNX(1,IE) DNX, DNY, DNZをDNXYZに融合 FX(1,IP2)=FX(1,IP2)-SWRK*DNXYZ(2,IE) FXYZ(1,IP2)=FXYZ(1,IP2)-SWRK*DNX(2,IE) FXYZ(2,IP1)=FXYZ(2,IP1)-SWRK*DNY(1,IE) FXYZ(2,IP1)=FXYZ(2,IP1)-SWRK*DNXYZ(2,1,IE) FX(1,IP2)=FX(1,IP2)-SWRK*DNX(2,IE) 1 FX(1,IP3)=FX(1,IP3)-SWRK*DNX(3,IE) FX(1,IP3)=FX(1,IP3)-SWRK*DNXYZ(3,IE) FXYZ(1,IP3)=FXYZ(1,IP3)-SWRK*DNX(3,IE) FXYZ(3,IP1)=FXYZ(3,IP1)-SWRK*DNZ(1,IE) FXYZ(3,IP1)=FXYZ(3,IP1)-SWRK*DNXYZ(3,1,IE) FX, FY, FZをFXYZに融合 2 FX(1,IP4)=FX(1,IP4)-SWRK*DNX(4,IE) FX(1,IP4)=FX(1,IP4)-SWRK*DNXYZ(4,IE) FXYZ(1,IP4)=FXYZ(1,IP4)-SWRK*DNX(4,IE) FXYZ(1,IP2)=FXYZ(1,IP2)-SWRK*DNX(2,IE) FXYZ(1,IP2)=FXYZ(1,IP2)-SWRK*DNXYZ(1,2,IE) 3 パターン2＋演算淳子変更 演算順序 パターンNo L1D ピーク キャッシュ 性能比% ミス率% FXYZ(2,IP2)=FXYZ(2,IP2)-SWRK*DNY(2,IE) FXYZ(2,IP2)=FXYZ(2,IP2)-SWRK*DNXYZ(2,2,IE) FY(2,IP1)=FY(2,IP1)-SWRK*DNXYZ(1,IE) FXYZ(2,IP1)=FXYZ(2,IP1)-SWRK*DNY(1,IE) FY(2,IP1)=FY(2,IP1)-SWRK*DNY(1,IE) パターン3＋パターン1 4 FXYZ(3,IP2)=FXYZ(3,IP2)-SWRK*DNZ(2,IE) FY(2,IP2)=FY(2,IP2)-SWRK*DNXYZ(2,IE) FXYZ(2,IP2)=FXYZ(2,IP2)-SWRK*DNY(2,IE) FXYZ(3,IP2)=FXYZ(3,IP2)-SWRK*DNXYZ(3,2,IE) FY(2,IP2)=FY(2,IP2)-SWRK*DNY(2,IE) FY(2,IP3)=FY(2,IP3)-SWRK*DNXYZ(3,IE) FXYZ(2,IP3)=FXYZ(2,IP3)-SWRK*DNY(3,IE) FY(2,IP3)=FY(2,IP3)-SWRK*DNY(3,IE) FXYZ(1,IP3)=FXYZ(1,IP3)-SWRK*DNX(3,IE) FXYZ(1,IP3)=FXYZ(1,IP3)-SWRK*DNXYZ(1,3,IE) FY(2,IP4)=FY(2,IP4)-SWRK*DNY(4,IE) FY(2,IP4)=FY(2,IP4)-SWRK*DNXYZ(4,IE) FXYZ(2,IP4)=FXYZ(2,IP4)-SWRK*DNY(4,IE) FXYZ(2,IP3)=FXYZ(2,IP3)-SWRK*DNY(3,IE) FXYZ(2,IP3)=FXYZ(2,IP3)-SWRK*DNXYZ(2,3,IE) 要求バイトが変わらないため、理論性能は FXYZ(3,IP3)=FXYZ(3,IP3)-SWRK*DNZ(3,IE) FZ(3,IP1)=FZ(3,IP1)-SWRK*DNZ(1,IE) FZ(3,IP1)=FZ(3,IP1)-SWRK*DNXYZ(1,IE) FXYZ(3,IP1)=FXYZ(3,IP1)-SWRK*DNZ(1,IE) FXYZ(3,IP3)=FXYZ(3,IP3)-SWRK*DNXYZ(3,3,IE) FZ(3,IP2)=FZ(3,IP2)-SWRK*DNXYZ(2,IE) FXYZ(3,IP2)=FXYZ(3,IP2)-SWRK*DNZ(2,IE) FZ(3,IP2)=FZ(3,IP2)-SWRK*DNZ(2,IE) 変化せず FXYZ(1,IP4)=FXYZ(1,IP4)-SWRK*DNX(4,IE) FXYZ(1,IP4)=FXYZ(1,IP4)-SWRK*DNXYZ(1,4,IE) FZ(3,IP3)=FZ(3,IP3)-SWRK*DNXYZ(3,IE) FXYZ(3,IP3)=FXYZ(3,IP3)-SWRK*DNZ(3,IE) FZ(3,IP3)=FZ(3,IP3)-SWRK*DNZ(3,IE) FXYZ(2,IP4)=FXYZ(2,IP4)-SWRK*DNY(4,IE) FXYZ(2,IP4)=FXYZ(2,IP4)-SWRK*DNXYZ(2,4,IE) FZ(3,IP4)=FZ(3,IP4)-SWRK*DNXYZ(4,IE) FXYZ(3,IP4)=FXYZ(3,IP4)-SWRK*DNZ(4,IE) FZ(3,IP4)=FZ(3,IP4)-SWRK*DNZ(4,IE) FXYZ(3,IP4)=FXYZ(3,IP4)-SWRK*DNZ(4,IE) FXYZ(3,IP4)=FXYZ(3,IP4)-SWRK*DNXYZ(3,4,IE) 1 2 34 45 3.95 4.05 4.05 4.41 3.98 3.69 3.69 3.37 メモリ スルー プット GB/s 35.70 35.91 35.88 38.80 37

スレッド並列化の⽅法の⼯夫
• 矩形分割ではブロック間に要素数のイン
バランスが生ずる可能性
• 要素のブロック化にMETISを用いて様々
な形状にも対応

オリジナル
形状

均等な要素
を含む複数
ブロックに
分割

)*+,-. /

("!!#
'"!!#
&"!!#
%"!!#

$"!!#
!"!!#
$!!!

各ブロックを
カラーリング
ブロック内には
たくさんの要素

ピーク性能⽐4.58%
L1Dミス率3.41%
メモリスループット
41.22GB/s

ピーク性能⽐
$!!!!

$!!!!!

$!!!!!!

012345678+ 9:4;<=

38

ཧ↷ࡍࡿࡼ࠺ 1 ḟඖ┠ࡀඛ࡟ືࡃࡼ࠺₇⟬㡰ᗎࢆ 㸬ࡍ࡞ࢃࡕ࣮࢝ࢿࣝࡢせồ B/F ್ࡣኚ໬ࡏࡎ⌮ㄽ まとめ ᑡࡋࡓࡓࡵ㸪࢟ࣕࢵࢩࣗࣛ࢖ࣥࡢ➇ྜࡀపῶࡉࢀࡓࡓࡵࣃ ࢱ࣮ࣥ 2/35.0ࡼࡾࡶࡉࡽ࡟㧗࠸ᛶ⬟ࢆ♧ࡋࡓ࡜⪃࠼ࡽࢀࡿ㸬 5.0 83%࡛࠶ࡿ㸬 4.0 ᅗ 18 ࡟ࡇࡇࡲ࡛࡛᭱Ⰻࡢ⤖ᯝࢆᚓࡓ㓄ิ⼥ྜࣃࢱ࣮ࣥ 4 命令実行 4.0 YZ(1,IP1)=FXYZ(1,IP1)-SWRK*DNX(1,IE) 待& YZ(2,IP1)=FXYZ(2,IP1)-SWRK*DNY(1,IE) YZ(3,IP1)=FXYZ(3,IP1)-SWRK*DNZ(1,IE) 3.0 YZ(1,IP2)=FXYZ(1,IP2)-SWRK*DNX(2,IE) YZ(2,IP2)=FXYZ(2,IP2)-SWRK*DNY(2,IE) 2.0 YZ(3,IP2)=FXYZ(3,IP2)-SWRK*DNZ(2,IE) Time[s] 3.0 ࡢྛࢫࣞࢵࢻࡢᐇຠ᫬㛫ෆヂࢆ♧ࡍ㸬ᅗ 2࿨௧ᐇ⾜ ࡟♧ࡋࡓ࢜ࣜࢪ '()*+,-関連待& Time [s] ࠾ࡾ㸪࢟ࣕࢵࢩ࣓ࣗࣔࣜࡢ฼⏝ຠ⋡ࡀᝏ࠸ࡇ࡜ࡀ࣎ࢺࣝࢿ ෆ࣮ࣝࣉ IE ࡢᅇ㌿ෆ࡛࢔ࢡࢭࢫࡉࢀࡿ㓄ิᩘࡀ኱ࡁࡃῶ ࢵࢡ࡟࡞ࡗ࡚࠸ࡿ࡜᥎ ࡉࢀࡿ㸬 䝞䝸䜰ྠᮇᚅ䛱 ࢼࣝࢥ࣮ࢻ࡜ẚࡋ࡚ィ⟬᫬㛫඲య࡟༨ࡵࡿ࣓ࣔࣜ࢔ࢡࢭࢫ 2.0 㛵㐃ࡢᚅࡕ࡜ࣂࣜ࢔ྠᮇᚅࡕࡢ๭ྜࡀ኱ᖜ࡟๐ῶࡉࢀ㸪࡯ 䝯䝰䝸䜰䜽䝉䝇 1.0 㛵㐃ᚅ䛱 ࡜ࢇ࡝ࡢ᫬㛫ࢆ࿨௧ᐇ⾜ࡀ༨ࡵࡿࡼ࠺࡟࡞ࡗࡓ㸬඲యࡢᐇ 0.0 YZ(1,IP3)=FXYZ(1,IP3)-SWRK*DNX(3,IE) 1.0 YZ(2,IP3)=FXYZ(2,IP3)-SWRK*DNY(3,IE) YZ(3,IP3)=FXYZ(3,IP3)-SWRK*DNZ(3,IE) 1 2 3 4 5 6 7 8 ຠ᫬㛫ࡶ 4.72 ⛊࠿ࡽ 1.76 ⛊࡬࡜ 2.7 ಸྥୖࡋࡓ㸬 䝇䝺䝑䝗 0.0 ケース名 内容 ᅗ 17 ࣃࢱ࣮ࣥ 3 㓄ิ FX, FY, FZ ࡢ⼥ྜղ ORG オリジナル (1,IP1)=FXYZ(1,IP1)-SWRK*DNXYZ(1,1,IE) 節点リナンバー，要素ブロック， A (2,IP1)=FXYZ(2,IP1)-SWRK*DNXYZ(2,1,IE) 外側カラー化 (3,IP1)=FXYZ(3,IP1)-SWRK*DNXYZ(3,1,IE) B 配列融合 (1,IP2)=FXYZ(1,IP2)-SWRK*DNXYZ(1,2,IE) (2,IP2)=FXYZ(2,IP2)-SWRK*DNXYZ(2,2,IE) C ブロック化手法改善 (3,IP2)=FXYZ(3,IP2)-SWRK*DNXYZ(3,2,IE) 2.0 C Sec.9 ᅗ 5 ࢫࣞࢵࢻẖࡢᐇ⾜᫬㛫ෆヂ 5. ࣓ࣔࣜ࢔ࢡࢭࢫ㣕ࡧࡢపῶ 1.5 ࿨௧ᐇ⾜ Time[s] YZ(1,IP4)=FXYZ(1,IP4)-SWRK*DNX(4,IE) ORG A B ORG Sec.7 Sec.8 YZ(2,IP4)=FXYZ(2,IP4)-SWRK*DNY(4,IE) YZ(3,IP4)=FXYZ(3,IP4)-SWRK*DNZ(4,IE) ࢟ࣕࢵࢩ࣓ࣗࣔࣜࡢ฼⏝ຠ⋡ࡀᝏ࠸ཎᅉ࡜ࡋ࡚⪃࠼ࡽ 1.0 䝞䝸䜰ྠᮇᚅ䛱 ࢀࡿࡢࡣ㸪᏶඲࡟⌮᝿ⓗ࡞≧ែ࡛ࡣ඲࡚ L1 ࢟ࣕࢵࢩࣗୖ 0.5 䝯䝰䝸䜰䜽䝉䝇 ࡟஌ࡿ࡜ぢ࡞ࡏࡿࡓࡵ↓ど࡛ࡁࡿࡣࡎࡢ㓄ิ FX, FY, FZ ࡬ 㛵㐃ᚅ䛱 0.0 ࡢ࢔ࢡࢭࢫࡀ㸪ࡇࡇ࡛⏝࠸࡚࠸ࡿ⌧ᐇⓗ࡞ධຊࢹ࣮ࢱ࡛ࡣ 1 2 3 4 5 6 7 8 䝇䝺䝑䝗 ⌮᝿≧ែ࠿ࡽ኱ࡁࡃࡎࢀ࡚࠾ࡾ㸪ࡑࡢࡓࡵ࣓ࣔࣜ࢔ࢡࢭࢫ 39

流体解析ソルバー ストアのシーケンシャル化

チューニング対象ルーチン • メインループ（時間積分ループ）の内訳 33% • 入力データは単純な6面体要素の形状で，1プロセスあたり 100万要素 • 全体では24プロセスで，京で1ノード1プロセス，1プロセス8ス レッド実行 • GRAD, DIV) 質量保存（ポアソン方程式）関係 • MOM, AXmult) 運動量関係 • CONV) データの変換 • COMM) 通信 • MISC) 色々雑多 • OTHER, 個別に扱えない細かいものの合計 一番重たいGRADを対象とする 41

チューニング対象カーネル Initialize array to zero Element Wise Product of vector (Hadamard Product) A={a1, a2, …, an} B={b1, b2, …, bn} C={c1, c2, …, cn} ci = ai bi Data Structure Conversion Main Calculation Ingredient of routine GRAD 42

カーネル処理内容 1 Elem.1 1 Vertex この節点は 要素１の頂点３ 要素２の頂点１ 要素３の頂点２ 要素４の頂点１ Pressure P1 rtex Ve 3 Elem.4 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 rtex Ve 1 DO iCOL=1,iNUM_COLOR 2 DO iBLK=1,iNUM_BLOCK(iCOL) Multi-Thread 3 iHEAD=iHEAD_ELEM(iBLK,iCOL) 4 iTAIL=iTAIL_ELEM(iBLK,iCOL) 5 DO iE=iHEAD, iTIAL 6 iN1=iNODE(1,iE) 7 iN2=iNODE(2,iE) ・・・ ・・・ 13 iN8=iNODE(8,iE) 14 fP=fPRES(iE) 15 16 fDP(1,iN1)+=fP*fDN(1,1,iE) 17 fDP(2,iN1)+=fP*fDN(2,1,iE) 18 fDP(3,iN1)+=fP*fDN(3,1,iE) 19 20 fDP(1,iN2)+=fP*fDN(1,2,iE) 21 fDP(2,iN2)+=fP*fDN(2,2,iE) 22 fDP(3,iN2)+=fP*fDN(3,2,iE) ・・・ ・・・ 44 fDP(1,iN8)+=fP*fDN(1,8,iE) 45 fDP(2,iN8)+=fP*fDN(2,8,iE) 46 fDP(3,iN8)+=fP*fDN(3,8,iE) 47 ENDDO 48 ENDDO 49 ENDDO 2 2 Vertex 前ページの事例の最終状態(のすこし違うVer) x2 Verte Vertex 3 4 Ni,j は有限要素法の形状関数で，要素iの頂点jごと定義されている， たとえば，N1,1 と N4,3 は一見同じ位置にあるが，違うもの 中央の節点における圧力勾配は以下の計算となる 43 43

データ構造 Elements Region Blocked Colored 1 5 2 6 9 13 10 14 3 7 4 8 11 15 12 16 Elements were sorted and renumbered in data Red Elements 1 2 3 Blue Yellow … 100 101 102 103 … 200 … 401 402 403 … 500 … 801 802 803 … 900 … 120112021203 … 1300 … Blk.1 iHEAD for Blk.1 Green Blk.2 iTAIL for Blk.1 Blk.5 Blk.9 Blk.13 44

性能評価 Runnig Time Performance Ratio to 128GFLOPS Memory throughput 91.5 sec 6.6 GFLOPS 5.2 % 26.8 GB/S Peak value is 46GB/s by STREAM bench. 前ページの事例と同じように，B/Fから理想的な性能を求め ると京では13% (16.7GFLOPS/NODE)出るとの見積 見積方法) 最内ループ1回転，つまり1要素当たり132Byte使う．一方演算量は48FLOPなので B/F=132/48=2.75 実効のメモリバンド幅はストリームベンチマークから46GB．演算性能は128GFLOPS/ノードなの で，B/Fの上限はB/F=46/128=0.36 その比から 0.36/2.75=13% 45

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 46

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 47

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 48

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 49

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 50

何がボトルネックか︖ 1 Elem.1 1 Vertex ¶N ì ¶N ü ÑP1 + = í 1,1 P1 , 1,1 P1 ý ¶y î ¶x þ Vertex 1 2 Vertex x2 Verte ①For Elem.1 2 Vertex 2 ¶N ì ¶N ü ÑP2 + = í 1, 2 P1 , 1, 2 P1 ý ¶y î ¶x þ Vertex 3 ¶N ì ¶N ü ÑP3 + = í 1,3 P1 , 1,3 P1 ý ¶y î ¶x þ rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 ②For Elem.2 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 プロファイラで実行時間内訳をみてみると • メモリアクセス待ち時間が多かった • ストア命令がたくさん出ていた ¶N 2,1 ü ì ¶N ¶N 2, 2 ü ì ¶N ÑP3 + = í 2,1 P2 , P2 ý ÑP2 + = í 2, 2 P2 , P2 ý ¶ x ¶ y ¶y î þ î ¶x þ ¶N 2,3 ü ì ¶N ÑP4 + = í 2,3 P2 , P2 ý ¶y î ¶x þ ③For Elem.3 ・・・略・・・ ④For Elem.4 ストアが多かった理由 ¶N 4,1 ü ì ¶N ÑP3 + = í 4,1 P4 , P4 ý ¶y î ¶x þ • ひとつの節点を複数の要素で共有していて，違う タイミングで同じ節点に何度もストアされる ¶N 4,3 ü ì ¶N ÑP1 + = í 4,3 P4 , P4 ý ¶y î ¶x þ ¶N 4, 2 ü ì ¶N ÑP5 + = í 4, 2 P4 , P4 ý ¶y î ¶x þ 節点への足しこみ(A=A+B)があるので • SIMDとソフトウェアパイプラインが出来ない 51

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて節点に代入(A=B) 52

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて要素に代入(A=B) 53

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて要素に代入(A=B) 54

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて要素に代入(A=B) 55

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて要素に代入(A=B) 56

ストアを連続，かつ1回だけにする 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte rtex Ve Pressure P1 3 Elem.4 rtex Ve 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 Vertex 3 4 処理の流れを変更する 要素を順番に巡って，要素からの計算結果を， 要素が使っている節点に足しこむ(A=A+B) ↓ 節点を順番に巡って，その節点を使っている 要素からの計算結果をまとめて要素に代入(A=B) 57

ソース変更 ループ構造を変更した 要素→節点と辿るのではなく節点のループを回して直接節点を参照する 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte Pressure P1 rtex Ve 3 Elem.4 2 5 Elem.2 Elem3 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 rtex Ve 1 DO iN=1,iNUM_NODE Multi-Thread 2 fTMPX=0.0 3 fTMPY=0.0 4 fTMPZ=0.0 5 DO iI=1,8 節点iNを使っている上位要素の数だけ回る 6 iE=iELEM(iI,iN) 上位要素番号 7 iV=iERT(iI,iN) 要素iEから⾒た節点iNは頂点番号 8 fP=fPRES(iE) 9 10 fTMPX += fP * fDN(1, iV, iE) 11 fTMPY += fP * fDN(2, iV, iE) 12 fTMPZ += fP * fDN(3, iV, iE) 13 ENDDO 14 fDP(1,iN) = fTMPX 15 fDP(2,iN) = fTMPY 16 fDP(3,iN) = fTMPZ 17 ENDDO Vertex 3 4 ケース1) 処理の流れだけ変更 58

ソース変更 もう1ケース ループ構造を変更した 要素→節点と辿るのではなく節点のループを回して直接節点を参照する 1 Elem.1 1 Vertex 2 2 Vertex x2 Verte Pressure P1 rtex Ve 3 Elem.4 2 5 Elem.2 Pressure P2 Elem3 Vertex 1 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 rtex Ve 1 DO iN=1,iNUM_NODE Multi-Thread 2 fTMPX=0.0 3 fTMPY=0.0 4 fTMPZ=0.0 5 DO iI=1,8 6 iE=iELEM(iI,iN) 7 fP=fPRES(iE) 8 9 fTMPX += fP * fDN2(1, iI, iN) 10 fTMPY += fP * fDN2(2, iI, iN) 11 fTMPZ += fP * fDN2(3, iI, iN) 12 ENDDO 13 fDP(1,iN) = fTMPX 14 fDP(2,iN) = fTMPY 15 fDP(3,iN) = fTMPZ 16 ENDDO Vertex 3 4 などの値は節点が持つようにした ケース2) 配列fDN→fDN2と変えた 59

改善効果 測定結果 1 A 123 5.6 4.4 23.2 2 B 135 5.2 4.1 23.2 A 44.3 15.8 12.4 47.2 単位 B 44.2 秒 15.9 GFLOPS 12.4 % 47.2 GB/S ケース2はほぼ理論的な性能12%に到達．一方でケース1は低い(メモリスループットも) 人工的なリスト値(全て同じ値)での測定結果 ケース パターン 実行時間 性能 ピーク比 メモリスループット 1 A 37.4 18.9 14.8 30.9 2 B 30.2 23.3 18.2 38 A 40.9 17.1 13.4 48.4 91秒 上 2倍向 リスト値 iE=iELEM(iI, iN) iV=iVERT(iI, iN) iI, iNの値によらず同じ値とし た ↓ キャッシュミスしなくなり，遅い 原因がキャッシュミスなのか否 かが明らかになる ケース パターン 実行時間 性能 ピーク比 メモリスループット 単位 B 41.7 秒 16.9 GFLOPS 13.2 % 47.7 GB/S 44秒 リスト値(配列iELEM,iVERT)の値を全て同一値にしてリストの性質の善し悪しの影響を排除した結果， 1の低性能はリストの影響と判明．2は元々リストアクセスされる配列が少ないので，リスト値の品質の 影響は小さい ※K.Kumahata, K.Minami, Y.Yamade, C.Kato, “Performance improvement of the general-purpose CFD code Front/Flow/blue on the K computer”, HPC Asia 2018, Tokyo, Japan, (2018) https://doi.org/10.1145/3149457.3149470 60

同じ改善を別カーネルへ アプリ中に同じ改善が適応可能なカーネルがあり，それにも適⽤し，実⾏時間半減 CASE ASIS tune1A tune2B tune2B1 tune2C 向上率 実行時間 1以上なら チューニング内容 [秒] 高速化 32.0 1.00 オリジナル(要素で回すループ) 55.1 0.58 節点で回すループに変更（配列アクセスはランダム) 19.1 1.67 tune1Aを配列を並び替えて配列アクセスを連続アクセス化 20.7 1.54 tune2Bで最内ループ回転数固定化 14.0 2.29 tune2B1で最内ループ内処理のアンロール実装を再度ループ化 ソース概要 60.0 Time [s] 50.0 40.0 30.0 20.0 10.0 0.0 ASIS tune1A tune2B tune2B1 tune2C DO IE=1, 全ヘキサ要素 DO I=1, ヘキサ要素の8頂点 IP=頂点Iのグローバル節点番号 演算(長いので略) 節点IPへのストア ENDDO ENDDO 配列へのストアは連続に、かつA=A+Bの形にしない 61

流体解析ソルバー ストリーム削減・セクタキャッシュ・プリフェッチ

チューニング対象カーネル FLD3X2ソース 1 Elem.1 1 Vertex Pressure P1 3 Elem.4 Elem.2 rtex Ve Elem3 2 Vertex 1 Pressure P2 Pressure P3 x3 Verte Vertex 3 Vertex Vertex 3 節点5 1 Pressure P4 1 Vertex 2 5 ※) 絵では3頂点の要素だが実際は8頂点の要素 2 2 Vertex x2 Verte rtex Ve 1 DO iE=1,iNUM_ELEM Multi-Thread 2 fLP(iE)=fDNX(1,iE)*fDP(1,iNODE(1,iE))+ 3 fDNX(2,iE)*fDP(1,iNODE(2,iE))+ ・・・ ・・・ 9 fDNX(8,iE)*fDP(1,iNODE(8,iE))+ 10 11 fDNY(1,iE)*fDP(2,iNODE(1,iE))+ 12 fDNY(2,iE)*fDP(2,iNODE(2,iE))+ ・・・ ・・・ 18 fDNY(8,iE)*fDP(2,iNODE(8,iE))+ 19 20 fDNZ(1,iE)*fDP(3,iNODE(1,iE))+ 21 fDNZ(2,iE)*fDP(3,iNODE(2,iE))+ ・・・ ・・・ 27 fDNZ(8,iE)*fDP(3,iNODE(8,iE)) 28 ENDDO 処理内容 要素が頂点として参照する節点がもつ(前の事例 GRAD求められた)圧力勾配ベクトルの発散を求 める Vertex 3 4 63

チューニング対象カーネル FLD3X2ソース 1 DO iE=1,iNUM_ELEM Multi-Thread 2 fLP(iE)=fDNX(1,iE)*fDP(1,iNODE(1,iE))+ 3 fDNX(2,iE)*fDP(1,iNODE(2,iE))+ ・・・ ・・・ 9 fDNX(8,iE)*fDP(1,iNODE(8,iE))+ 10 11 fDNY(1,iE)*fDP(2,iNODE(1,iE))+ 12 fDNY(2,iE)*fDP(2,iNODE(2,iE))+ ・・・ ・・・ 18 fDNY(8,iE)*fDP(2,iNODE(8,iE))+ 19 20 fDNZ(1,iE)*fDP(3,iNODE(1,iE))+ 21 fDNZ(2,iE)*fDP(3,iNODE(2,iE))+ ・・・ ・・・ 27 fDNZ(8,iE)*fDP(3,iNODE(8,iE)) 28 ENDDO 処理内容 要素が頂点として参照する節点がもつ(前の事例 GRAD求められた)圧力勾配ベクトルの発散を求 める 要素の頂点数 fLP fDNX fDP fDNZ fDP 2 1 fLP(1) fLP(4) 節点5 5 ※) 絵では3頂点の要素だが実際は8頂点の要素 fDNY fDP fLP(2) 3 fLP(3) 4 64

性能評価 FLD3X2ソース 1 DO iE=1,iNUM_ELEM Multi-Thread 2 fLP(iE)=fDNX(1,iE)*fDP(1,iNODE(1,iE))+ 3 fDNX(2,iE)*fDP(1,iNODE(2,iE))+ ・・・ ・・・ 9 fDNX(8,iE)*fDP(1,iNODE(8,iE))+ 10 11 fDNY(1,iE)*fDP(2,iNODE(1,iE))+ 12 fDNY(2,iE)*fDP(2,iNODE(2,iE))+ ・・・ ・・・ 18 fDNY(8,iE)*fDP(2,iNODE(8,iE))+ 19 20 fDNZ(1,iE)*fDP(3,iNODE(1,iE))+ 21 fDNZ(2,iE)*fDP(3,iNODE(2,iE))+ ・・・ ・・・ 27 fDNZ(8,iE)*fDP(3,iNODE(8,iE)) 28 ENDDO 最適化状況 ストアの向きがGRADとは逆なのでデータ依存性 なしでSIMD化，ソフトウェアパイプライン化は適応 されている 結果 実行時間 性能 ピーク比 メモリスループット 57.4 秒 10.3 GFLOPS 8.0 % 36.8 GB/S B/Fから求まる理想値は13%なので遅い 65

改善⼿法1 FLD3X2ソース 1 DO iE=1,iNUM_ELEM Multi-Thread 2 fLP(iE)=fDNX(1,iE)*fDP(1,iNODE(1,iE))+ 3 fDNX(2,iE)*fDP(1,iNODE(2,iE))+ ・・・ ・・・ 9 fDNX(8,iE)*fDP(1,iNODE(8,iE))+ 10 11 fDNY(1,iE)*fDP(2,iNODE(1,iE))+ 12 fDNY(2,iE)*fDP(2,iNODE(2,iE))+ ・・・ ・・・ 18 fDNY(8,iE)*fDP(2,iNODE(8,iE))+ 19 20 fDNZ(1,iE)*fDP(3,iNODE(1,iE))+ 21 fDNZ(2,iE)*fDP(3,iNODE(2,iE))+ ・・・ ・・・ 27 fDNZ(8,iE)*fDP(3,iNODE(8,iE)) 28 ENDDO 改善手法 1．配列のサイズをコンパイラに教える 上位ルーチンでは real*4 :: fDNX(8、NE) call FLD3X2(・・・,fDNX,8,NE,・・・) FLD3X2ルーチンでは subroutine FLD3X2(・・・,fDNX,SZ,NE,・・・) real*4 :: fDNX(SZ,NE) となっており、コンパイラは配列fDNXの1次元⽬の ⼤きさが可変という前提で命令を作るが、実際には 固定なので、コンパイラに教えてみる 66

改善⼿法2 FLD3X2ソース 1 DO iE=1,iNUM_ELEM Multi-Thread 2 fLP(iE)=fDNX(1,iE)*fDP(1,iNODE(1,iE))+ 3 fDNX(2,iE)*fDP(1,iNODE(2,iE))+ ・・・ ・・・ 9 fDNX(8,iE)*fDP(1,iNODE(8,iE))+ 10 11 fDNY(1,iE)*fDP(2,iNODE(1,iE))+ 12 fDNY(2,iE)*fDP(2,iNODE(2,iE))+ ・・・ ・・・ 18 fDNY(8,iE)*fDP(2,iNODE(8,iE))+ 19 20 fDNZ(1,iE)*fDP(3,iNODE(1,iE))+ 21 fDNZ(2,iE)*fDP(3,iNODE(2,iE))+ ・・・ ・・・ 27 fDNZ(8,iE)*fDP(3,iNODE(8,iE)) 28 ENDDO 改善手法 2. SIMDロードを適応するため，アンロール実装 をループに戻す 1 iE+1 iE fDNX 2 3 4 5 6 7 8 1 2 ひとつのロード命令でまとめてロードできるのだが、 fDNX(1,iE)のような書き方だとやってくれないので、 fDNX(j, iE)のように書けるように変更(リループ) 67

改善⼿法2 DO iE=1,iNUM_ELEM スレッド並列、SWP DO J=1,8 SIMD, FULLUNROLL tmp = tmp + fDNX(J,iE)* fDP(1, iNODE(J,iE)) + fDNY(J,iE)* fDP(2, iNODE(J,iE)) + fDNZ(J,iE)* fDP(3, iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO 配列の1次元目がループ変数になったため、SIMDロードになるはず 68

改善⼿法3 FLD3X2ソース 改善手法 DO iE=1,iNUM_ELEM スレッド並列、SWP 3. DO J=1,8 SIMD, FULLUNROLL tmp = tmp + fDNXYZ(J,1,iE)* fDP(1, iNODE(J,iE)) + fDNXYZ(J,2,iE)* fDP(2, iNODE(J,iE)) + fDNXYZ(J,3,iE)* fDP(3, iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO do I=i, N A(i) = A(i) + B(i) * C(i) enddo Register/演算器 A(i) ＝ αB(i) ＋ C(i) do I=i, N A(i) = A(i) + D(1,i) * D(2,i) enddo B(i) C(i) Aは書き込みだけでなく 参照もされている場合 Cache 演算結果 Memory Register/演算器 cache line ③ ロード ④ ① B cache line ② C cache line 配列fDNX, fDNY, fDNZをfDNXYZにまとめ てストリーム削減 A 1 1 2 Cache 演算結果 A(i) ＝ αB(i) ＋ C(i) D B C B C B C B C B C B ・・・ 2 3 3 4 4 5 5 6 Memory cache line Bn Cn Bm Cm ② ロード ③ ① A Bn Cn Bm Cm 69

改善⼿法4 FLD3X2ソース 改善手法 !ocl cache_sector_size(3,9) !ocl cache_subsector_assign(fDP) DO iE=1,iNUM_ELEM スレッド並列、SWP DO J=1,8 SIMD, FULLUNROLL tmp = tmp + fDNXYZ(J,iE)* fDP(1, iNODE(J,iE)) + fDNXYZ(J,iE)* fDP(2, iNODE(J,iE)) + fDNXYZ(J,iE)* fDP(3, iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO 4. セクタキャッシュを利用 Way 0 Line 1 Line 2 Line N 1 2 3 4 5 6 7 8 先 優 9 10 11 P D f ・・・・・・・・・・・・・・・・・・・・・・・・ 配列fDPは繰り返し利用されると分かっているの で優先的にキャッシュに乗っていて欲しい L2は12ウェイあるので、3：9に分けて、9の側は fDPを優先にする 2 fDP(1:3,2) fDP(1:3,1) 1 fLP(1) fLP(4) 節点5 fLP(2) 3 fDP(1:3,3) 5 fDP(1:3,5) fLP(3) fDP(1:3,4) 4 70

改善効果 オリジナル 実行時間 性能 ピーク比 メモリスループット 改善手法 57.4 秒 10.3 GFLOPS 8.0 % 36.8 GB/S 1. 配列のサイズをコンパイラに教える 2. 1＋SIMDロードを適応するため，アンロール実装 をループにする 3. 2＋配列fDNX，fDNY，fDNZをfDNXYZにまとめて ストリーム削減 4. 3＋配列を「京」の機能であるセクタキャッシュに載 せる 1 実行時間 性能 ピーク比 メモリスループット 2 3 4 52.6 49.9 49 47.5 秒 11.3 12.4 12.6 13.0 GFLOPS 8.8 9.7 9.8 10.2 % 40.3 42.4 43.2 42.5 GB/S 71

ソフトウェアプリフェッチ＠富岳試作機 !ocl cache_sector_size(3,9) !ocl cache_subsector_assign(fDP) DO iE=1,iNUM_ELEM スレッド並列、SWP DO J=1,8 SIMD, FULLUNROLL tmp = tmp + fDNXYZ(J,1,iE)* fDP(1,iNODE(J,iE)) + fDNXYZ(J,2,iE)* fDP(2,iNODE(J,iE)) + fDNXYZ(J,3,iE)* fDP(3,iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO 富岳試作機で測定 演算量は全てのスレッドで⼀定のはずが、バリア同期が発⽣ している ↓ ⼊⼒メッシュデータに依存して、スレッドごとのキャッシュミス率が 異なってしまっている 72

ソフトウェアプリフェッチ＠A64FX !ocl cache_sector_size(3,9) !ocl cache_subsector_assign(fDP) DO iE=1,iNUM_ELEM スレッド並列、SWP オリジナル !ocl prefetch_read(fDP(1,NODE(1,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(2,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(3,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(4,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(5,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(6,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(7,iE+32)),level=1) !ocl prefetch_read(fDP(1,NODE(8,iE+32)),level=1) DO J=1,8 SIMD, FULLUNROLL tmp = tmp + fDNXYZ(J,1,iE)* fDP(1,iNODE(J,iE)) + fDNXYZ(J,2,iE)* fDP(2,iNODE(J,iE)) + fDNXYZ(J,3,iE)* fDP(3,iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO 配列fDPのプリフェッチをパラメータスタディ スレッドごとの打ち分けが均等になり、1秒縮んだ 73

命令 SIMD化 4命令 時 減少 ソフトウェアプリフェッチ＠A64FX „ !ocl cache_sector_size(3,9) !ocl cache_subsector_assign(fDP) DO iE=1,iNUM_ELEM スレッド並列、SWP 総命令数 削減 オリジナル DO J=1,8 SIMD, FULLUNROLL !ocl prefetch_read(fDP(1,NODE(J,IE+64), level=1) tmp = tmp + fDNXYZ(J,1,iE)* fDP(1,iNODE(J,iE)) + fDNXYZ(J,2,iE)* fDP(2,iNODE(J,iE)) + fDNXYZ(J,3,iE)* fDP(3,iNODE(J,iE)) ENDDO fLP(iE) = tmp ENDDO プリフェッチの指⽰⾏をSIMDが有効になっている ループ内へと変更し、プリフェッチもSIMD化 10.1秒 → 9.2秒 → 6.7秒 74

乱流燃焼解析コード

チューニング対象計算機 FX1000 Clock Freq. 2.2 GHz Core GFLOPS 70 GFLOPS L1D Cache 64 KiB (256GB /s) CMG CPU # of Core 12 GFLOPS 884 GFLOPS Memory 8 GiB (256 GB/s) L2 Cache 8 Mib (1024 GB/s) # of CMG 4 GFLOPS 3379 GFLOPS Memory 32 GiB Node # of CPU # of Node Syste GFLOPS m Memory 1 リングバスには他にPCIeやインターコネクトのコントローラが接続 1コア当たりの演算性能 2.2GHz x 2FMA x 8SIMD × 2Pipe = 70.4GFLOPS Machine 京 5760 SIMD width 2 (128 bit) 19.4 PFLOPS # of register 180 TiB 256 FX100 FX1000 4 (256 bit) 8 (512 bit) 256 32 • SIMDありきの性能なのに注意 • レジスタ数が減っているのに注意 76

ベンチマーク問題 液体酸素と⽔素の同軸噴流⽕炎 反応機構で素反応を含む8化学種 H2 150K Cost distribution on SORA by profiler FIPP Use 1CMG (12threads) on FX1000 100step 10MPa O2 100K 7100 Cells (1cell for thickness) p=2 (# of SPs in Cell 27) 実在気体の状態⽅程式 SRK⽅程式 輸送係数 燃焼モデル A Update T and P H Correct Residual Chung B Compute E I Compute Residual F.G. CHEMKIN C Compute Cp J Compute Self Resiaual(A) Wilke混合則 D Compute Trans Prop K Compute Self Residual(B) E Compute Cv L Compute Gradient F Compute Thermal Prop M Sum of under 1% G Ref. flamelet tab. Flamelet Progress Variable (FPV) Procedures for Combustion Calculations 77

インライン展開のやりすぎに注意 チューニング) ⼀般的には演算区間であるような下位のループの実⾏時間 を短縮し，GFLOPSやGB/sを向上させること． D) Compute Trans Propのコールツリー 通らない(実⾏時 パラメータ依存) インライン展開 本コードでは，上位ループから回転数・演算数がごく⼩さい 下位ループを多数呼んでおり，⼿続き呼び出しのオーバー ヘッドや，レジスタ不⾜に伴う演算待ちがボトルネックであり， 現状，それらの低減が主眼． ⾮インライン展開 ⼿法ごとに整理して述べる デフォルト 最良 これ以外のループの回転数はごく少ない (例:コンパイル時に明らかになっている数8) do icell=1, セル数 スレッド並列化 do isp=1, 27 Call_Subroutines enddo 最も回転数が多いループicell enddo 2番⽬に回転数が多いループisp 原点が0でないので注意 ComputeMwBar ComputeRconst FX100向けにインライン展開を多⽤．使⽤するレジスタ数 が増えるが，FX1000ではレジスタ数が減っているので，レ ジスタに乗せるデータの⼀時退避・復帰が頻発し性能劣化 ComputeTransProp_LowPress ComputeTransProp_Chung_Single Mixing_Wilke 0ならインライン展開する，1ならインライン展開しない 78

⼩回転ループはさけたい • コンパイル時にはすでに決まっている，⼩回転数のループ • ソフトウェアパイプライン化されるが，回転数が少ないため効果なし • ループ制御のオーバーヘッドを避けるためループで無くした⽅がよい F) Compute Thermal Propの内の下位のループ SIMD化(VL: 8) do icell=1, セル数 スレッド並列化 SOFTWARE PIPELINING化 do isp=1, 27 IPC: 3.00, ITR: 144, 今のケースでは実際は2 ・・・ MVE: 5, POL: S) do ichem=1, pr%nchem y(ichem) = value(ie + ichem, icell) enddo ・・・ コンパイル時には数が分かっているが enddo 可変回転数のループとして扱われている enddo プリ処理部でpr%nchem=_NSCLと設定 このようなループがある FULL UNROLLING do icell=1, セル数 スレッド並列化 do isp=1, 27 ・・・ do ichem=1, _NSCL y(ichem) = value(ie + ichem, icell) enddo ・・・ enddo enddo 原点が0でないので注意 79

⼩回転ループはさけたい 別のルーチンにも同様の⼿法を当てはめてみた Loop B,C) ソース上で直接回転数をハードコー ドしている do i=1,10など 52.6 52.5 52.4 52.3 52.2 52.1 52.0 Loop C Time [s] Loop A) 前ページと同じく，コンパイル時には 回転数が確定しているが，コンパイ ラにはわからない 今のケースでは8回転 0) Original (Implemented as Loop) 1) Modified (Unrolled) 0 1 0 1 0 1 0 1 Loop B 0 0 1 1 0 0 1 1 Loop A 0 0 0 0 1 1 1 1 原点が0でないので注意 80

OpenMPのスケジュール調整 A) Update T and P ルーチンの実⾏時間内訳 Inctruction Commit Time [s] 10.0 8.0 Barrier Wait 6.0 4.0 do itr=1, max res = ・・・ if(res < EPS) exit enddo Instruction Wait 2.0 0 1 2 3 4 5 6 7 8 9 10 11 Thread Memory/Cach e Acces s W ait ■Number of FLOP 1.0E+11 2.0E+10 8.0E+10 1.5E+10 6.0E+10 1.0E+10 4.0E+10 2.0E+10 5.0E+09 0.0E+00 0.0E+00 0 1 2 3 4 5 6 7 8 9 10 11 Thread 収束判定に合格して ループを抜けるタイミングが まちまちだった • 当該ループのOpenMPのディレクティブに schedule(dynamic, 1) を指定し !$omp parallel do schedule(dynamic, 1) とした ■Number of Instructions 0.0 • FX1000のCPU性能解析レポートの出⼒，バリア同期待ちが多い • セルを巡るループでスレッド並列化しており，各スレッドへのセルの割当量 は均⼀だが演算数に違いが出ている 12.0 Inctruction Commit 10.0 Time [s] 12.0 8.0 Barrier Wait 6.0 4.0 Instruction Wait 2.0 0.0 0 1 2 3 4 5 6 7 8 9 10 11 Thread Memory/Cach e Acces s W ait 81

重複した処理の回避 A) Update T and Pルーチン中の温度更新反復のコールツリー • 同じサブルーチンが複数の上位から呼ばれていることが多い Update T and P • 例えば⾚枠ConvertMass2Moleは質量分率からモル分率に変 換するルーチン．反復ループ外から1回，反復ループ中からは他 のルーチンを介して1反復当たり3回呼ばれるが，全ての呼び出し が同じ計算 • 他にもループ中に温度に依存しない計算があったり，依存するが 同じ温度から同じ結果得る計算を何度も読んでいる個所があり， 除去を⾏った このとき浮動⼩数点演算数は 1.09×1012 から 8.10x1011 に25％ほど削 減 温度を更新する 反復ループがある do itr=1, iter_max call subroutine T = ... enddo 82

配列初期化時の性能劣化 ※前ページまでの事例とは，ソースの版数と⼊⼒データが異なる C ase N o. A T otal range1 当該ルーチン概要 01 subroutine sample(n1, n2, input_arrayA, input_arrayB) 02 implicit none 03 1ルーチン単体(下位のルーチン含まず) 04 ! argument でコストの18%弱を占める 05 integer, intent(in) :: n1, n2 06 real(8), intent(in) :: input_arrayA(n1), input_arrayB(n2) ルーチン内に細かくタイマールーチンを挿⼊． 07 特定部分を計るためのタイマーが有効な 08 ! local 09 real(8) :: work_arrayA(n1), work_arrayB(n2) 時だけ速い 10 ・ 11 ・ 12 ! call timer_start(range2) 13 work_arrayA(:) = input_arrayA(:) 14 ! call timer_stop(range2) 15 16 ! call timer_start(range3) 17 work_arrayB(:) = input_arrayB(:) range3 20 ! call timer_stop(range3) 21 ・ 22 ・ range2 range3 range4 other 23 end subroutine sample 31 12.80 0.12 32 2.81 0.11 0.12 33 2.99 0.11 0.11 34 2.81 0.11 1.08 11.60 0.11 2.48 • range2,3は引数の値をワーク配列へコピーする配列式 0.10 0.11 2.56 • 0.11 0.11 2.49 コンパイラは配列式をループとして扱いループ最適化 (実⾏時回転数889以上でスレッド並列，SIMD，ソフトウェアパイプライン など適応) • コメントアウトしてあるタイマールーチン呼び出しを有効にすると速くなる 単位 [s] range2,3でおかしなことが起きている 83

配列初期化時の性能劣化 コンパイル時及び実⾏時の以下の条件が重なり遅くなっていた 1. コンパイラは配列式をループとして扱うため，range2とrange3は それぞれループとなる 2. ⾃動スレッド並列化が有効である(-Kparallel)場合，同じ ループのスレッド並列版と，逐次版を⽣成し，実⾏時のループ回 転数で切り替える 3. スレッド⽣成のオーバーヘッドを避けるため，並列リージョン拡張と いう，可能なら複数のループを合成してからスレッド並列化の処 理を⾏う機能がある この機能が動作したときは，逐次版のループを作らない range2と3の間にタイマールーチン呼び出しがある場合は並列 リージョン拡張は働かない このような現象が出るときの解決策 • -Knoregion_extensionを指定 • 配列式に!OCL SERIALを指定 • 今回のように無害なルーチンを配列式の間に挿⼊ range2 イメージ do i=1, n1 work_arrayA(i) = input_arrayA(i) enddo range3 do i=1, n2 work_arrayB(i) = input_arrayB(i) enddo イメージ range2 if(n1 > 889) then !$omp parallel range2のスレッド並列版ループ !$omp end parallel else range2の逐次版ループ endif イメージ !$omp parallel range2のスレッド並列版ループ range3のスレッド並列版ループ !$omp end parallel range3 if(n2 > 889) then !$omp parallel range3のスレッド並列版ループ !$omp end parallel else range3の逐次版ループ endif 実⾏時の回転数は range2は72, range3は9． 並列リージョン拡張されなかった場合は，逐 次版ループを実⾏する 拡張された場合は低回転ループにも関わら ず，スレッド並列版ループを実⾏する そのためオーバーヘッドで性能劣化 84

まとめ 100.0 アプリ全体の実⾏時間は，81秒程度から51秒程度に 改善 80.0 配列初期化時の性能劣化は，統⼀して測定できな かったので，その分10秒を個別に加えると90秒程度か ら51秒程度に改善 Time [s] 今回ご紹介した紹介したチューニング事例を全て適応 With Array Problem WithInit. Array Init Problem 60.0 Barrier Wait 40.0 Instruction Wait 20.0 0.0 As-Is Current Thread プロファイラ(CPU性能解析レポート)を使うとバーヘッドが出る ため，ここではプロファイラを⽤いて求めた実⾏時間の内訳を， プロファイラなしで実⾏した際の全体実⾏時間で補正 他にも，コンパイルオプションの再考として，元々FX100向け だったので⾊々指定されていた -Kfast,parallel,noalias=s,array_private, preex,ocl,nounroll,autoobjstack を再整理し て -Kfast,parallel,ocl として56.0秒から53.8秒 Inctruction Commit Memory/Cache Acces s W ait チューニング） ⼀般的には演算区間であるような下位のループの実⾏時 間を短縮し，GFLOPSやGB/sを向上させること． 本コードでは，上位ループから回転数・演算数がごく⼩さい 下位ループを多数呼んでおり，⼿続き呼び出しのオーバー ヘッドや，レジスタ不⾜に伴う演算待ちがボトルネックであり， 今回はそれらボトルネックの低減が主眼でした 85

その他のTips

なるべく⻑いループでSIMDに SIMD化は最内ループに適応される (1,2) do idx=1, max(=191700) SWP (3-1) do ichem=1, spc_num(=8) 8SIMD (3-2) do j=1,3 FULL UNROLL (3-3) do j=1, spc_num(=8) 8SIMD 1.03秒，3.33E+10 FLOP, SIMD率40.0% (1,2) do idx=1, max(=191700) 8SIMD (3-1) do ichem=1, spc_num(=8) FULL UNROLL (3-2) do j=1,3 FULL UNROLL (3-3) do j=1, spc_num(=8) FULL UNROLL 0.46秒，2.35E+10FLOP，SIMD率90.4% 87

Z-Fillの効果 Z-Fillなし Register/演算器 Cache 4stream ① cache line A(i) ＝ αB(i) ＋ C(i) cache line cache line Memory 暗黙のロード A ④ ② B ③ C 演算結果 Z-Fillあり -Kzfill キャッシュライン確保命令 Register/演算器 Cache 3stream cache line A(i) ＝ αB(i) + C(i) cache line cache line ① ② ③ Memory A B C 演算結果 88

Z-Fillの効果 疑似コード integer*8 :: I, SZ, ITR, NITR parameter(SZ=104857600) parameter(NITR=100) real*8 :: A(SZ), B(SZ), C(SZ), alpha Z-Fillにより暗黙のAのロードがなくなり ストリームが削減し性能改善 配列ひと つのサイ start = gettime() call FAPP_START() do ITR=1, NITR do I=1, SZ A(I)=alpha * B(I) + C(I) enddo 本体部分 enddo call FAPP_STOP() dt = start – gettime() case GB dt [s] GB/s ズ MB 800 80 8 withoutzfill 335.544 1.616 207.657 withoutzfill_numactl 335.544 1.618 207.343 withzfill 251.658 1.236 203.634 withzfill_numactl 251.658 1.237 203.382 withoutzfill 335.544 1.601 209.577 withoutzfill_numactl 335.544 1.623 206.767 withzfill 251.658 1.242 202.633 withzfill_numactl 251.658 1.241 202.788 withoutzfill 335.544 1.463 229.422 withoutzfill_numactl 335.544 1.674 200.500 withzfill 251.658 1.279 196.833 withzfill_numactl 251.658 1.284 195.948 89

その他 ifを含むループはリストベクトル化 オリジナルは，全ての節点で 判定とループを実行 上位要素数8までの節点 上位要素数8以上の節点 節点リスト 1 2 3 4 5 6 7 8 7.8秒 IF(NEP(8)>8) IF(NEP(8)>8) IF(NEP(8)>8) IF(NEP(8)>8) THEN IF(NEP(8)>8) THEN IF(NEP(8)>8) THEN IF(NEP(8)>8) THEN IF(NEP(8)>8) THENTHENTHENTHEN DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) DO I=9,NEP(8) ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ・・・ ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO ENDDO ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF ENDIF リストベクトル化 ifに合格するものだけの集合 上位要素数が8以上の節点リスト 3 4 6 配列、特に多次元配列の初期化を、配列の 形に合わせたようなループでやらない 8 IP=LST(4) IP=LST(4) IP=LST(4) IP=LST(4) DO I=9,NEP(IP) DO I=9,NEP(IP) DO I=9,NEP(IP) DO I=9,NEP(IP) ・・・ ・・・ ・・・ ・・・ ENDDO ENDDO ENDDO ENDDO ENDIF ENDIF ENDIF ENDIF 0.0?秒 リストベクトル化後は IFなしで，実行可能 ↓ 最内ループは必ず実行 され，SIMD化・SWP化 されている REAL A(N1, N2, NE) DO IE=1, NE マルチスレッド DO I=1, N2 DO J=1, N1 SIMD A(J,I,IE)=Const ENDDO ENDDO ENDDO 2.2秒 C言語のmemsetを呼ぶラッパー関数を 用意 配列を8等分し，8スレッドから呼ぶ REAL A(N1, N2, NE) 0.6秒 SZ=(N1*N2*NE) / 8 DO I=8 マルチスレッド ST=1+SZ*(I-1) call memset_wrap(A(P),Const,SZ) ENDDO 90

まとめ • これまでに単体性能チューニングを⾏った流体コード2種についての事例を紹介しました • 流体解析ソルバーでは • • キャッシュ利⽤効率の改善・スレッド並列化の⽅法の⼯夫 • ストアのシーケンシャル化 • ストリーム削減・セクタキャッシュ・プリフェッチ 乱流燃焼コードでは • インライン展開のやり過ぎに注意 • ⼩回転ループはさけたい • OpenMPのスケジュール調整 • 重複した処理の回避 • 配列初期化時の性能劣化 • その他細かいTips • 皆さんのプログラムく⾼速化の参考になると幸いです 91