並列前処理手法と領域分割， マルチコア時代の戦略

1. 並列前処理手法と領域分割，マルチコア時代の戦略並列前処理手法と領域分割，マルチコア時代の戦略 分野横断型研究会「アルゴリズムによる計算科学の融合と発展」 2009年4月22・23日　筑波大学計算科学研究センター 中島 研吾 東京大学情報基盤センター 海洋研究開発機構地球シミュレータセンター 科学技術振興機構戦略的創造研究推進事業（CREST）

2. 科学技術計算の真髄：SMASH • 幅広い分野 • バランス • 分野間協力 • 少しは他のカテゴリーについても知らないと協力は進まない。 • 応用 • 自分の問題が解ければOK • 巷に役に立つライブラリが無い（特に疎行列） • アルゴリズム • 万能を目指したい

3. 並列有限要素法による大規模シミュレーション並列有限要素法による大規模シミュレーション • 非構造格子 • 要素単位のローカルな処理⇒大規模疎行列 • 悪条件問題（ill-conditioned problems） • 前処理付並列反復法

4. 講演の概要 • 並列反復法と領域分割：SMASH • 選択的オーバーラップ〔KN 2007〕 • 階層型領域間境界分割（Hierarchical Interface Decomposition）〔Henon & Saad 2007〕 • 拡張型HID法の提案：SMASH • 悪条件向け領域分割手法 • Hybrid並列プログラミングモデル：SMASH • HIDと並列多重格子法（時間があれば）：SMASH

5. 「並列」前処理手法の技術的課題（の一部） • 領域分割 • 簡単な問題は簡単に解ける，効率は出やすい • 難しい問題はやっぱり難しい • Block Jacobi型局所前処理 ⇒ 領域数増加による反復回数増加 • 領域外の影響を（基本的に）無視 • 悪条件問題で顕著

6. 本講演で扱う悪条件問題 • 様々な悪条件問題がある • 普通の工学的な問題は大抵「悪条件問題」 • 係数行列の固有値分布，条件数 • ここでは，特に三次元固体力学における問題を扱う • 接触条件 • 不均質性 • 捩れ等 • Block ILU型の前処理手法 • 各節点に3自由度（変位3成分）：効率，安定性

7. 悪条件問題の例Heterogeneous Fields, Distorted Meshes

8. 地震発生サイクルシミュレーションにおける接触問題地震発生サイクルシミュレーションにおける接触問題 • プレート境界における準静的応力蓄積過程 • 非線形接触問題をNewton-Raphson法によって解く • ALM法（Augmented Lagrangean,拡大ラグランジェ法）による拘束条件：ペナルティ数 • 領域分割による並列有限要素法

9. 領域間境界＝接触面収束は最悪

10. 「硬い」要素群上に領域境界が来ると収束は悪化する「硬い」要素群上に領域境界が来ると収束は悪化する E=100 E=103 3D Solid Mechanics E: Young’s Modulus

11. 不均質弾性問題：203要素■：4×4×20

12. 不均質弾性問題：203要素BILU(0)-GPBiCG，反復回数不均質弾性問題：203要素BILU(0)-GPBiCG，反復回数 • ■■：n= 0.25 • ■：E=1.00 • 1-プロセッサ • ■：E=10-3， 34回 • ■：E=100 ， 31回 • ■：E=10+3, 84回 • 8-プロセッサ（オーバーラップ無し） • ■：E=10-3， 53回（×1.56） • ■：E=100 ， 52回（×1.68） • ■：E=10+3，158回（×1.88）

13. 悪条件問題への対処 • 悪条件問題を前処理付き反復法で解く場合，並列計算時には収束性が著しく悪化する場合があり，安定した収束を与える前処理，領域分割の研究は重要 • 領域分割：いわゆる数値計算ライブラリではカバーされていない分野 • 対処法 • マルチレベル解析，Coarsegrid法 • 深い領域間オーバーラップ • 計算・通信コスト

14. 領域間オーバーラップの拡張 Cost for computation and communication may increase

15. Overlap 深さ E=10-3 E=100 E=10+3 0 53 52 158 1 34 33 103 2 32 32 100 3 30 32 97 4 31 31 82 1領域 34 31 84 不均質弾性問題：203要素ILU(0)-GPBiCG，8領域，反復回数オーバーラップ領域拡張の影響

16. 接触問題：オーバーラップ拡張の効果 • [KN 2005] • BILU(0,1,2) • for “consistent” node number cases • IBM SP3 in NERSC/LBNL

17. 最近やっていること • 安定で効率的な並列前処理手法，領域分割手法開発 • アプリケーションの特徴を最大限使用 • 安定化，効率：ブロック化 • アプリケーション⇒特殊前処理⇒一般化 • 問題に応じて，最適な前処理手法，領域分割，パラメータの組み合わせを自動的に選択するための手法の確立 • 大域的・局所的情報の利用 • 大域的情報：係数行列の条件数等 • 局所的情報：各要素のローカルな情報 • アプリケーションの性質 • 実問題の特性を反映させたベンチマーク

18. 実問題の特性を反映させたベンチマーク • 疎行列解法，前処理手法の検証 • Matrix Market，実アプリマトリクス • 制約が多い⇒結局ごく限られた条件を代表 • 取得が困難な場合がある • 長時間のシミュレーション，大規模マトリクスデータ • 非線形問題：違うフェーズ⇒性質違うマトリクス • 実問題の特性を反映させたベンチマーク • 元の問題と類似した性質の係数マトリクスを生成 • 元の問題が非線形でも，何らかの形で線形化が行われているはず • 「ベンチマーク」は線形問題でも良い • パラメータ，形状，問題規模を自由に変えられる • 係数マトリクスが導出された過程くらいはわかっていなくてはならない

19. 領域分割手法 • 選択的オーバーラップ 〔KN 2007〕 • 様々な条件に対応できるよう，有限要素法によるアプリケーションの特性（要素の属性，物性）を利用した前処理手法，領域分割手法 • 適応的にオーバーラップ深さを調整 • 選択的フィルイン：前処理手法 • 階層型領域間境界分割（HID）〔Henon & Saad 2007〕 • Hierarchical Interface Decomposition • 領域数増加による反復回数増加の効率的な抑止 • PHIDAL（Parallel Hierarchical Interface Decomposition Algorithm ）アルゴリズム

20. 選択的オーバーラップ，選択的フィルイン • アセンブリ構造物における接触問題 • アプリケーションの特性を利用 • 選択的フィルイン：前処理 • 選択的オーバーラップ：領域分割

21. アセンブリ構造の例：ジェットエンジン

22. 接触面における節点整合，不整合 接触面 整合 不整合 アセンブリ構造：部位ごとに別々に メッシュを作るのでこのようなことが おこりうる 地震：大すべり

23. 例題：接触面節点不整合問題アセンブリ構造を模擬例題：接触面節点不整合問題アセンブリ構造を模擬 • 各ブロックは1辺長さ1.0の立方体要素（ヤング率=1.00，ポアソン比=0.30の弾性体）に分割。 • 各ブロックは0.10ずつ離れており，その間が交差するトラス要素（弾性体）で結合されている。

24. 例題：接触面節点不整合問題アセンブリ構造を模擬例題：接触面節点不整合問題アセンブリ構造を模擬 • トラス要素のヤング率をブロック部分の103とすることによって接触面における拘束条件を模擬。 • z=0でz方向の変位を固定，x=0及びy=0で対称とし，z=zmaxの面にz方向に一様分布荷重を与えている。

25. Selective fill-ins • トラス要素に接続している節点のみに1レベル高いfill-inを適用する • 本問題の場合： BILU(1+) • Block ILU（三次元弾性問題，1節点3自由度） • BILU(2)：トラス要素に接続する節点 • BILU(1)：それ以外の節点 • 計算量はBILU(1)並みであるが，前処理性能としてはBILU(2)相当が期待される。

26. ● 2nd order fill-in’s are considered for these nodes ● 2nd order fill-in’s are NOT considered for these nodes ● 2nd order fill-in’s are NOT considered for these nodes Idea of “Selective fill-ins”: ILU(1+)

27. Selective Overlapping • 「Selective fill-ins」の考え方を領域間オーバーラップに拡張する • 一般の節点に関しては領域間オーバーラップの拡張を「遅らせる」 • 領域間オーバーラップ拡張による計算量，通信量の増加を抑制できる。

28. Internal Nodes for Partitioning ● Internal Nodes Domain Boundary

29. One-Layer Overlapping (d=0/1) ● Internal Nodes ● External Nodes ■ Overlapped Elements This is the general configuration of local data set for parallel FEM (one-layer of overlapping).

30. Extension of Overlapped Zones (2-layers: d=2) ● Internal Nodes ● External Nodes ■ Overlapped Elements

31. Extension of Overlapped Zones Extension of Overlapped Zones (d=2 and d=1+) ● Internal Nodes ● External Nodes ■ Overlapped Elements

32. Extension of Overlapped Zones Extension of Overlapped Zones (d=2 and d=1+) ● Internal Nodes ● External Nodes ■ Overlapped Elements Selective Overlapping (d=1+) “Delayed” extension for elements which do not include nodes connected to truss-type elements

33. Extension of Overlapped Zones Extension of Overlapped Zones (d=3 and d=2+) ● Internal Nodes ● External Nodes ■ Overlapped Elements Selective Overlapping (d=2+) Reduced cost for computations and communications delayed delayed

34. BILU with selective fill-in/overlapping • BILU (p)-(d) • p level of fill-ins (0, 1, 1+, 2, 2+ …) • d depth of overlapping (0, 1, 1+, 2, 2+ …)

35. Results: 64 cores3,090,903DOF, l=103, e=10-8Effect of Overlapping ITERATIONS ● BILU(1)-(d) ■ BILU(1+)-(d) ▲ BILU(2)-(d) Elapsed Time #-non-zero’s of [M]

36. Results: 64 cores3,090,903DOF, l=103, e=10-8Effect of Overlapping ● BILU(1)-(d) ■ BILU(1+)-(d) ▲ BILU(2)-(d) ITERATIONS Elapsed Time

37. HID：Hierarchical Interface Decomposition〔Henon & Saad〕 • 階層型領域間境界分割 • 二言くらいで言うと： • 階層的な領域分割 • Nested Dissection • 各レベルでは各領域は直接結合しない⇒レベル内並列性 • 計算コスト的には（d=0）と（d=1）の中間くらい：低コスト • d：オーバーラップ

38. Parallel ILU for each Connector at each LEVEL • レベルの若い順に番号を振りなおす • 各レベル内で不完全LU（0）分解を並列に実施可能 • sub-domainの番号を0,1,2,3とする • 数字は隣接するsub-domainの番号

39. 各レベルにおける通信の発生 • 「高い」レベルのコネクタに属する節点群における計算 • 「低い」レベルのコネクタの計算結果を利用 • 計算済 • 隣接する「低い」レベルのコネクタに属する節点群が他の領域に属している場合 • 通信が必要となる

40. Forward Substitution 余計な通信 が発生

41. 計算結果（64 cores）接触問題3,090,903DOF ■BILU(p)-(0): Block Jacobi ■BILU(p)-(1) ■BILU(p)-(1+) ■BILU(p)-HID GPBiCG

42. 0 1 2 3 Level-1 0,1 0,2 2,3 1,3 0,1, 2,3 Level-2 Level-4 HID vs. Selective Overlapping • HIDと選択的オーバーラップはほぼ同じ性能だが後者が若干良い • 特に条件の悪い問題 • BILUで高いorderのfill-inが必要となるような問題 • オリジナルのHIDはfill-inのorderが高くなると，領域外の（同じレベルにある節点の）fill-inの影響を考慮できない • ILU(0)としては完全だが・・・

43. HID vs. Selective Overlapping • HIDと選択的オーバーラップはほぼ同じ性能だが後者が若干良い • 特に条件の悪い問題 • BILUで高いorderのfill-inが必要となるような問題 • オリジナルのHIDはfill-inのorderが高くなると，領域外の（同じレベルにある節点の）fill-inの影響を考慮できない • ILU(0)としては完全だが・・・ 0 1 Level-1 2 3 0,1 0,2 Level-2 2,3 1,3 0,1, 2,3 Level-4

44. 並列反復法と領域分割 • 選択的オーバーラップ〔KN 2007〕 • 階層型領域間境界分割（Hierarchical Interface Decomposition）〔Henon & Saad 2007〕 • 拡張型HID法の提案 • 悪条件向け領域分割手法 • Hybrid並列プログラミングモデル • HIDと並列多重格子法（時間があれば）

45. 要素がねじれた問題（1/3） • 3D linear elastic problem with locally distorted elements

46. 要素がねじれた問題（2/3） • 3D linear elastic problem with locally distorted elements • 立方体メッシュ • Z軸周りに回転 • 局所的な不均質性 • 局所的なねじれの程度を表す • sequential Gauss algorithm [Deutsch & Journel 1988]

47. 要素がねじれた問題（3/3） • 3D linear elastic problem with locally distorted elements

48. ■BILU(1)-(d,a) GPBiCG ■BILU(1+,120°)-(d,a) ■BILU(1+, 60°)-(d,a) ■BILU(1+, 30°)-(d,a) ■BILU(2)-(d,a) 要素がねじれた問題（不均質分布）BILU(p,q)-(d,a)3,090,903DOF，64コアMAX distortion: 150-deg. (q,a)

49. ■BILU(1)-(d,a) GPBiCG ■BILU(1+,120°)-(d,a) ■BILU(1+, 60°)-(d,a) ■BILU(1+, 30°)-(d,a) ■BILU(2)-(d,a) 要素がねじれた問題（不均質分布）BILU(p,q)-(d,a)3,090,903DOF，64コアMAX distortion: 225-deg. (q,a)

50. HIDの改良を試みるExtended Version of HID • オーバーラップ領域の拡張 • レベル間セパレータを「厚く」する • Thicker Separators