マルチコア時代の並列前処理手法 Parallel l Preconditioning i Methods for Iterative Solvers in Multi-Core Era 中島研吾 東京大学情報基盤センター 2010 年 10 月 18 日 京都大学数理解析研究所 (RIMS) 研究集会 : 科学技術計算アルゴリズムの数理的基盤と展開
2 We are now in Post Peta Scale Era PFLOPS: Peta (=10 15 ) Floating OPerations per Sec. Exa FLOPS (=10 18 ) will be attained in 2018 or 2019
3 Exa-Scale Systems Peta-scale -> Evolution, Exa-scale -> Revolution 様々な技術的問題点 ( 例 ) >10 8 コア数を持つシステムの耐故障性 (Fault Tolerance) 電力消費量 現状の最も効率的なシステム :2MW/PFLOPS( 年 2 億円 ) ExaFLOPS:2GW, 年 2,000 億円 20MW にすることが必要 メモリーウォール問題 現状 Byte/Flop rate (B/F) 0.40 -> 0.10, 0.02? 汎用的システムは困難 分野間協力重要 H/W, S/W, Applications 計算機科学, 計算科学, 数値アルゴリズム
4 IESP: International Exascale Software Project http://www.exascale.org/ International Project A single country cannot do that 4 Workshops since 2009 5 th is during October 18 th -19 th in Maui, HI, USA Current Status Discussions on Road-map
5 Key-Issues towards Appl./Algorithms on Exa-Scale Systems Jack Dongarra (ORNL/U. Tennessee) at SIAM/PP10 ( 日本応用数理学会誌 Vol.20-3に関連記事 ) Hybrid/Heterogeneous Architecture Multi + GPU Multi + Many (more intelligent) Mixed Precision Computation Auto-Tuning/Self-Adapting Adapting Fault Tolerant Communication Reducing Algorithms
ACES2010 Heterogeneous Architecture by 6 (CPU+GPU) or (CPU+Many) will be general in less than 5 years NVIDIA Fermi Intel Knights Ferry
ACES2010 7 CPU+Accelerator (GPU, Many) 高いメモリーバンド幅 現状の GPU には様々な問題点 通信 :CPU-GPU/GPU-GPU プログラミングの困難さ :CUDA,OpenCL O CL は状況を変えつつあるが 限定されたアプリケーションのみで高効率 : 陽的 FDM,BEM メニーコア (Manys) Intel Many Integrated Core Architecture (MIC) GPU より賢い : 軽い OS, コンパイラが使える Intel Knights Ferry with 32 s is available soon for use on development of programming environment (very limited users) Knights Corner with >50 s (22nm) in 2012 or 2013? 近い将来 GPU と Many(MIC 的な意味での ) は大差なくなる
8 Hybrid 並列プログラミングモデルは 必須 Message Passing MPI Multi Threading OpenMP
9 2010RIMS Flat MPI vs. Hybrid Flat-MPI:Each PE -> Independent memor ry memor ry memor ry Hybrid:Hierarchal Structure mem mory mem mory mem mory
2010RIMS 10 背景 T2Kオープンスパコンン ( 東大 ) 並列多重格子法 (Multigrid) 前処理付き CG 法 MGCG Flat MPI vs. Hybrid (OpenMP+MPI) Hybrid MPI のプロセス数を減らせる 通信オーバーヘッド減少 メモリ的には厳しくなる : 特に疎行列ソルバー
11 2010RIMS T2Kオープンスパコン仕様ン仕様 T2K( 東大 )(1/2) http://www.open-supercomputer.org/ 筑波大, 東大, 京大 T2Kオープンスパコン ( 東大 ) Hitachi HA8000クラスタシステム 2008 年 6 月 ~ 952ノード (15,232コア), 141 TFLOPS peak Quad- Opteron (Barcelona) TOP500 53 位 (Jun 2010)
12 2010RIMS T2K( 東大 )(2/2) AMD Quad- Opteron Memory Memory (Barcelona) 2.3GHz 4 sockets per node L2 L2 L2 L2 L1 L1 L1 L1 16 s/node L3 L3 L2 L2 L2 L2 L1 L1 L1 L1 マルチコア, マルチソケット cc-numa(cache coherent Non-Uniform Memory Core Core Core Core Core Core Core Core Access) L1 L1 L1 L1 L1 L1 L1 L1 ローカルメモリ上のデータをできるだけ使用する 陽的なコマンドラインスイッチ NUMA control Core Core Core Core L2 L2 L2 L2 L3 Core Core Core Core L2 L2 L2 L2 L3 Memory Memory
2010RIMS 13 Multigrid is scalable Weak Scaling: Problem Size/Core Fixed 三次元ポアソン方程式 ( 一様 ) 3000 2500 ICCG MGCG 2000 Iterations 1500 1000 500 0 1E+06 1.E+06 1E+07 1.E+07 1E+08 1.E+08 DOF
2010RIMS 14 Multigrid is scalable Weak Scaling: Problem Size/Core Fixed MGCG 法の計算時間は Weak Scaling では一定 =Scalable 3000 2500 ICCG MGCG 2000 Iterations 1500 1000 500 16 32 64 128 0 1E+06 1.E+06 1E+07 1.E+07 1E+08 1.E+08 DOF
15 2010RIMS Flat MPI vs. Hybrid Flat-MPI:Each PE -> Independent memor ry memor ry memor ry Hybrid:Hierarchal Structure mem mory mem mory mem mory
16 2010RIMS Flat MPI vs. Hybrid 性能は様々なパラメータの組み合わせによって決まる ハードウェア コア,CPUのアーキテクチュア ピーク性能 メモリ性能 ( バンド幅, レイテンシ ) 通信性能 ( バンド幅, レイテンシ ) それらのバランス アプリケーション 特性 :memory bound,communication bound 問題サイズ
2010RIMS 17 Flat MPI, Hybrid (4x4, 8x2, 16x1) Higher Performance of HB16x1 is important Flat MPI 0 1 2 3 Hybrid 4x4 0 1 2 3 Hybrid 0 1 2 3 8x2 Hybrid 0 1 2 3 16x1
2010RIMS 18 Domain Decomposition Inter Domain: MPI-Block Jacobi Intra Domain: OpenMP-Threads (re-ordering) example: 6 nodes, 24 sockets, 96 s Flat MPI HB 4x4 HB 16x1
2010RIMS 19 解析対象 透水係数が空間的に分布する三次元地下水流れ ポアソン方程式 透水係数は地質統計学的手法によって決定 Deutsch & Journel, 1998 規則正しい立方体ボクセルメッシュを使用した有限体積法 局所細分化を考慮 周期的な不均質性 : 128 3 φ φ λ + λ + x x y y φ = 0@ x = x max z φ λ z = q
Groundwater Flow through Heterogeneous Porous Media Homogeneous Uniform Flow Field Heterogeneous Random Flow Field
2010RIMS 21 前処理付き CG 法 Multigrid id 前処理 線形ソルバーの概要 IC(0) for Smoothing Operator (Smoother) Additive Schwartz Domain Decomposition 並列 ( 幾何学的 ) 多重格子法 当方的な8 分木 V-cycle 領域分割型 :Block-Jacobi 局所前処理, 階層型領域間通信 最も粗い格子 ( 格子数 =プロセッサ数 ) は1コアで実施
2010RIMS 22 IC(0) as smoother of Multigrid IC(0) is generally more robust than GS. IC(0) smoother with Additive Schwartz Domain Decomposition (ASDD) provides robust convergence and scalable performance of parallel computation, even for ill- conditioned problems KN 2002.
23 Overlapped Additive Schwartz Domain 2010RIMS pp Decomposition Method for Stabilizing Localized Preconditioning for Stabilizing Localized Preconditioning Global Operation Global Operation Ω Mz = r Local Operation Ω 1 Ω 2 2 2 2 1 1 1, Ω Ω Ω Ω Ω Ω = = r z M r z M n n Global Nesting Correction 2 2 2 1 1 1 Ω Ω Ω Ω Ω Ω Global Nesting Correction Ω 1 Ω 2 ( ) 1 1 1 1 1 2 1 2 1 1 1 1 1 1 Γ Γ Ω Ω Ω Ω Ω Ω + n n n n z M z M r M z z Γ 2 1 Γ 1 2 ( ) 1 1 1 1 2 1 2 1 2 2 2 2 2 2 Γ Γ Ω Ω Ω Ω Ω Ω + n n n n z M z M r M z z
2010RIMS T2K/Tokyo Hardware/Software up to 512 nodes (8,192 s) Program Hitachi FORTRAN90 + MPI CRS matrix storage CM-RCM Reordering for OpenMP Ax-b / b =10-12 for Convergence 不均質性 最大最小透水係数の比 = 10 10 (10-5 ~10 +5 ) Multigrid id Cycles 1 V-cycle/iteration for (i=0; i<n; i++) { for (k=index(i-1); k<index(i); k++{ Y[i]= Y[i] + A [k]*x[item[k]]; } } 2 smoothing iterations for restriction/prolongation at every level 1 ASDD iteration cycle for each resrtiction/prolongation 24
Algorithm09 25 前処理付き反復法の SMP/Multi での OpenMP による並列化 DAXPY, SMVP, Dot Products 簡単 前処理 :ILU 系分解, 前進後退代入 大域的な依存性 (Global dependency) 並び替え (Reordering) による並列性の抽出 Multicolor Ordering (MC), Reverse-Cuthill-Mckee (RCM) 同じ色内の要素は独立 並列化可能 地球シミュレータ 向け最適化 [KN 2002,2003] 並列及びベクトル性能 並列性高く安定な CM-RCM を採用
2010RIMS 26 Ordering Methods Elements in same color are independent: to be parallelized 29 45 22 61 16 46 11 62 47 7 63 4 48 2 64 1 29 22 16 11 7 4 2 1 53 36 19 2 49 33 17 1 37 49 13 30 29 23 51 14 17 30 12 53 15 31 8 55 16 5 32 3 37 30 23 17 12 8 5 3 7 54 37 20 3 50 34 18 44 41 38 57 31 42 24 58 18 43 13 59 44 9 60 6 44 38 31 24 18 13 9 6 25 8 55 38 21 4 51 35 50 33 9 45 25 39 35 10 32 26 25 37 11 19 27 14 39 12 10 28 50 45 39 32 25 19 14 10 43 26 9 56 39 22 5 52 55 37 51 53 46 38 40 54 33 39 26 55 20 40 15 56 55 51 46 40 33 26 20 15 61 44 27 10 57 40 23 6 59 17 5 56 21 52 19 6 47 22 41 21 7 34 23 27 23 8 21 24 59 56 52 47 41 34 27 21 14 62 45 28 11 58 41 24 62 33 60 49 57 34 53 50 48 35 42 51 35 36 28 52 62 60 57 53 48 42 35 28 31 15 63 46 29 12 59 42 64 1 63 17 61 32 58 18 54 53 49 19 43 74 36 20 64 63 61 58 54 49 43 36 48 32 16 64 47 30 13 60 MC (Color#=4) Multicoloring RCM Reverse Cuthill-Mckee CM-RCM (Color#=4) Cyclic MC + RCM
Effect of Optimization 64 s (4 nodes) of T2K/Tokyo 64 3 cells/ 16,777,216 cells Full Optimization NUMA Control First Touch Data Placement Further Reordering (with Contiguous/Sequential Memory Access)
Algorithm09 28 Policy ID Command line switches 0 no command line switches 1 2 3 4 --cpunodebind=$socket --interleave=all --cpunodebind=$socket --interleave=$socket --cpunodebind=$socket --membind=$socket --cpunodebind=$socket --localalloc 5 --localalloc l ll sec. 100.00 80.0 60.0 40.0 NUMA control Memory Memory L3 L3 L2 L2 L2 L2 L2 L2 L2 L2 L1 L1 L1 L1 L1 L1 L1 L1 Core Core Core Core Core Core Core Core Core Core Core Core Core Core Core Core L1 L1 L1 L1 L1 L1 L1 L1 L2 L2 L2 L2 L2 L2 L2 L2 L3 L3 Memory Memory Initial NUMA control Full Optimization Down is good 20.0 0.0 Flat MPI HB 4x4 HB 8x2 HB 16x1
Algorithm09 29 First Touch Data Placement 配列のメモリ ページ : 最初にtouchしたコアのローカルメモリ上に確保計算と同じ順番で初期化 do lev= 1, LEVELtot do ic= 1, COLORtot(lev)!$omp parallel l do private(ip,i,j,isl,iel,isu,ieu) i i i i i do ip= 1, PEsmpTOT do i = STACKmc(ip,ic-1,lev)+1, STACKmc(ip,ic,lev) RHS(i)= 0.d0; X(i)= 0.d0; D(i)= 0.d0 isl= indexl(i-1)+1 iel= indexl(i) do j= isl, iel iteml(j)= 0; AL(j)= 0.d0 enddo isu= indexu(i-1)+1 ieu= indexu(i) do j= isu, ieu itemu(j)= 0; AU(j)= 0.d0 enddo enddo enddo!$omp omp end parallel do enddo enddo
Further Re-Ordering for Continuous Memory Access: Sequential 5 colors, 8 threads Initial Vector Coloring (5 colors) +Ordering color=1 color=2 color=3 color=4 color=5 Coalesced (Original) i color=1 color=2 color=3 color=4 color=5 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 Sequential 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 8 8 8 8 8
2010RIMS 31 Flat MPI, Hybrid (4x4, 8x2, 16x1) Higher Performance of HB16x1 is important Flat MPI 0 1 2 3 Hybrid 4x4 0 1 2 3 Hybrid 0 1 2 3 8x2 Hybrid 0 1 2 3 16x1
Effect of F.T. + Sequential Data Access 16,777,216= 64x64 3 cells, 64 s, CM-RCM(2) Time for Linear Solvers 100.0 Initial NUMA control Full Optimization 80.0 se ec. 60.0 40.0 Down is good 20.0 0.0 Flat MPI HB 4x4 HB 8x2 HB 16x1
Effect of F.T. + Sequential Data Access 128 3 tri linear hexahedral elements, 6,291,456 DOF ICCG Solvers for 3D Linear Elastic Eqn s, 32 nodes of T2K (512 s), Time for Linear Solvers, HB 4x4 is the fastest UP is good 33 Rela ative Perf formance 1.50 1.00 0.50 000 0.00 Initial CASE-1 CASE-2 CASE-3 Flat MPI HB 4x4 HB 8x2 HB 16x1 coalesced coalesced + NUMA coalesced + NUMA+ first touch sequential + NUMA + first touch Parallel Programming Models
Effect of Number of Colors
色数の効果 (CM-RCM) 16,777,216= 64x64 3 cells, 64 s 色数が増えると収束は改善, 計算時間は CM-RCM(2) が最も短い Iterations sec. Iterations 70 28.0 65 60 55 50 45 Flat MPI HB 4x4 HB 8x2 HB 16x1 1 10 100 1000 COLOR# sec. 26.0 24.0 22.0 20.0 T2K: Flat MPI T2K: HB 4x4 T2K: HB 8x2 T2K: HB 16x1 1 10 100 1000 COLOR#
色数の効果 (CM-RCM) 16,777,216= 64x64 3 cells, 64 s 色数が増えると収束は改善, 計算時間は CM-RCM(2) が最も短い : 反復あたり計算時間短い sec./iter sec. 0.60 28.0 sec./iteratio on 0.50 040 0.40 0.30 T2K: Flat MPI T2K: HB 4x4 T2K: HB 8x2 T2K: HB 16x1 1 10 100 1000 COLOR# sec. 26.0 24.0 22.0 20.0 T2K: Flat MPI T2K: HB 4x4 T2K: HB 8x2 T2K: HB 16x1 1 10 100 1000 COLOR#
色数の効果 (CM-RCM) RCM: 前進後退代入時に変数値が変わるため, キャッシュラインからメモリに戻されてしまう可能性がある 29 22 16 11 7 4 2 1 37 30 23 17 12 8 5 3 44 38 31 24 18 13 9 6 50 45 39 32 25 19 14 10 55 51 46 40 33 26 20 15 59 56 52 47 41 34 27 21 62 60 57 53 48 42 35 28 45 10 39 5 35 2 33 1 17 46 11 40 6 36 3 34 53 18 47 12 41 7 37 4 24 54 19 48 13 42 8 38 59 25 55 20 49 14 43 9 29 60 26 56 21 50 15 44 63 30 61 27 57 22 51 16 61 29 62 30 63 31 64 32 25 57 26 58 27 59 28 60 53 21 54 22 55 23 56 24 17 49 18 50 19 51 20 52 45 13 46 14 47 15 48 16 9 41 10 42 11 43 12 44 37 5 38 6 39 7 40 8 64 63 61 58 54 49 43 36 32 64 31 62 28 58 23 52 1 33 2 34 3 35 4 36 RCM CM-RCM(2) MC(2)
2010RIMS 38 Weak Scaling Up to 8,192 s (512 nodes) 64 3 cells/ 2,147,483,648 cells CM-RCM(2)
2010RIMS 39 Weak Scaling 64 3 cells/, up to 8,192 s (2.05 10 9 cells) sec. Iterations sec. 100 80 60 40 Flat MPI init. HB 4x4 init. HB 8x2 init. HB 16x1 init. tions Itera 200 150 100 Flat MPI init. HB 4x4 init. HB 8x2 init. HB 16x1 init. 20 50 0 10 100 1000 10000 CORE# 0 10 100 1000 10000 CORE#
2010RIMS 40 Coarse Grid Solver の改良 領域数が増えると反復回数が増加 ( 特に Flat MPI) 最も粗い格子 (Coarse Grid Solver) Iteratio ons 200 150 100 各領域 1メッシュになった状態で1コアに集める 50 IC(0) スムージングを一回施す Coarse Grid Solver 改良 IC(0) スムージングを収束 (ε=10-12 ) まで繰り返す :C1 マルチグリッド (V-cycle) を適用し, 収束 (ε=10-12 ) まで繰り返す (8,192=32 16 16): C2 0 Flat MPI init. HB 4x44 init. it HB 8x2 init. HB 16x1 init. 10 100 1000 10000 CORE#
2010RIMS 41 Weak Scaling: Flat MPI 64 3 cells/, up to 8,192 s (2.05 10 9 cells) sec. Iterations 100 80 Flat MPI init. Flat MPI C1 Flat MPI C2 200 150 Flat MPI init. Flat MPI C1 Flat MPI C2 sec. 60 40 tions Itera 100 20 50 0 10 100 1000 10000 0 10 100 1000 10000 CORE# CORE#
2010RIMS 42 Weak Scaling: Flat MPI 64 3 cells/, up to 8,192 s (2.05 10 9 cells) Coarse Grid Solver Iterations grid solve er) sec c. (coarse 1.E+02 1E+01 1.E+01 1.E+00 1.E-01 1.E-02 Flat MPI init. Flat MPI C1 Flat MPI C2 tions Itera 200 150 100 50 Flat MPI init. Flat MPI C1 Flat MPI C2 1.E-03 10 100 1000 10000 CORE# 0 10 100 1000 10000 CORE#
2010RIMS 43 Weak Scaling: Flat MPI 64 3 cells/, up to 8,192 s (2.05 10 9 cells) sec. Iterations 40 35 Flat MPI init. Flat MPI C1 Flat MPI C2 100 90 Flat MPI init. Flat MPI C1 Flat MPI C2 sec. 30 tions 80 25 Iterat 70 20 60 15 10 100 1000 10000 CORE# 50 10 100 1000 10000 CORE#
2010RIMS 44 Weak Scaling 64 3 cells/, up to 8,192 s (2.05 10 9 cells) at 8,192 s: Flat MPI(35.7sec), HB 4x4(28.4), 8x2(32.8), 16x1(34.4) sec. Iterations sec. 100 80 60 40 Flat MPI C2 HB 4x4 C2 HB 8x2 C2 HB 16x1 C2 tions Itera 200 150 100 Flat MPI C2 HB 4x4 C2 HB 8x2 C2 HB 16x1 C2 20 50 0 10 100 1000 10000 CORE# 0 10 100 1000 10000 CORE#
2010RIMS 45 Strong Scaling 512x256x256= 33,554,432 cells Up to 1,024 s (64 nodes) CM-RCM(2)
Strong Scale: Parallel Performance 512x256x256= 33,554,432 cells based on performance of Flat MPI with 16 s HB 4x4 at 1,024 s: 73.7% Up is good Par rallel Perf formance (%) 120 100 80 60 40 20 Flat MPI HB 8x2 HB 4x4 HB 16x1 0 16 32 64 128 256 512 1024 CORE#
2010RIMS 47 関連研究 OpenMP/MPI Hybrid を並列多重格子法に適用した例は近年特に増加している : Sandia, LLNL Alison Baker (LLNL) et al., On the Performance of an Algebraic Multigrid Solver on Multi Clusters, (VECPAR 2010) Hypre Library (BoomerAMG), weak scaling Hera Cluster(T2K 東大とほぼ同じアーキテクチャ ) ~216 nodes, 3,456 コア ( 発表では >10,000 コア ) MultiCore SUPport library (MCSup) HB 4 4 が最も性能が良い
2010RIMS 48 まとめ ( 多重格子法 (MG) 前処理 +CG 法 ) 不均質多孔質媒体中の三次元地下水流れ, 有限体積法 IC(0) smoother + ASDD, 幾何学的 MG OpenMP/MPI Hybrid 並列プログラミングモデル on T2K ( 東大 ) NUMA Policy First Touch Data Placement + Sequential Reordering Coarse Grid Solver 改良 HB 4x4(a single MPI process per socket) が最も効率が良い : メモリを最も効率よく使っている, 通信オーバーヘッドも少ない 反復回数は並列プログラミングモデルによってほとんど変化しない Memory L3 L2 L2 L2 L2 L1 L1 L1 L1 Core Core Core Core Core Core Core Core L1 L1 L1 L1 L2 L2 L2 L2 L3 Memory Memory L3 L2 L2 L2 L2 L1 L1 L1 L1 Core Core Core Core Core Core Core Core L1 L1 L1 L1 L2 L2 L2 L2 L3 Memory
2010RIMS 49 今後の課題 粗い格子レベルにおけるコア数の漸減 全体のコア数, 領域数が増えると通信オーバーヘッドが増加 Hybrid における領域内並べ替え CM-RCM HID 並列化 : 結構時間がかかる Communication Reducing Algorithms 並列 MG: とにかく通信多い
Further Re-Ordering for Continuous Memory Access: Sequential 5 colors, 8 threads Initial Vector Coloring (5 colors)+ordering color=1 color=2 color=3 color=4 color=5 Coalesced (Original) color=1 color=2 color=3 color=4 color=5 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 1 2 3 4 5 6 7 8 Sequential 1 1 1 1 1 2 2 2 2 2 3 3 3 3 3 4 4 4 4 4 5 5 5 5 5 6 6 6 6 6 7 7 7 7 7 8 8 8 8 8