計算科学技術 特論B 2016年5月12日 国立研究開発法人 理化学研究所 計算科学研究機構�運用技術部門 ソフトウェア技術チーム�チームヘッド 南�一生 [email protected] RIKENADVANCEDINSTITUTEFORCOMPUTATIONALSCIENCE 2016年5月12日 計算科学技術 特論B 2016年5月12日 計算科学技術 特論B 2016年5月12日 計算科学技術 特論B 2016年5月12日 計算科学技術 特論B 5 2016年5月12日 計算科学技術 特論B A Processor Performance difficult High easy Low RSDFT Lattice QCD Seism3D NICAM PHASE FFB Low 2016年5月12日 計算科学技術 特論B 2016年5月12日 計算科学技術 特論B ParallelPerformance Performance Parallel Parallelization 7 difficult High ● ● 漏れ電流が問題 FET LSI FET Si n-Si LSI 2020 2030 LSI FET 漏れ電流を押さえる 2016年5月12日 計算科学技術 特論B 9 FET Kohn-Sham方程式 電子密度 n(r )=| ∑ ϕi (r ) 2 波動関数 i φi: 電子軌道(=波動関数) ハミルトニアン i:電子準位(=エネルギーバンド) r:空間離散点(=空間格子) Hϕ i (r) = ε iϕ i (r) 2016年5月12日 計算科学技術 特論B 10 実空間法 Hϕ i (r) = ε iϕ i (r) Kohn-Sham u 3 u ML2 ML1,ML2,ML3 ML(=ML1×ML2×ML3 3 2016年5月12日 計算科学技術 特論B 11 Self-ConsistentFieldprocedure 1 (CG)共役勾配法 2 (GS)Gram-Schmidt規格直交化 3 密度とポテンシャルの更新 4 (SD)部分対角化 SCF 2016年5月12日 計算科学技術 特論B 12 r u × × SCF 4 or 5 × SCF v Gram-Schmidt w n SCF n r o scf_mix um ※ → SCF r F → 2016年5月12日 計算科学技術 特論B u v END 13 J.-I.Iwataetal.,J.Comp.Phys.(2010) Realspace CPUspace Blue:Siatom Yellow:electrondensity 2016年5月12日 計算科学技術 特論B 14 2016年5月12日 計算科学技術 特論B RSDFT • 実空間差分法 • 空間並列 計算コアの最適化 • 行列積化 ターゲット計算機:PACS-CS,T2K-Tsukuba スレッド並列の実装 ターゲット計算機:PACS-CS,T2K-Tsukuba 2016年5月12日 計算科学技術 特論B 16 M EOT UP ψ 1" 4 ψ1 ψ 2" 4 ψ2 ψ 1" ψ 2 ψ 1" ψ 3" 4 ψ3 ψ 1" ψ 3 ψ 1" ψ 2" ψ 3 ψ 2" ψ 4" 4 ψ4 ψ 1" ψ 4 ψ 1" ψ 2" ψ 4 ψ 2" ψ 3" ψ 4 ψ 3" ψ 5" 4 ψ5 ψ 1" ψ 5 ψ 1" ψ 2" ψ 5 ψ 2" ψ 3" ψ 5 ψ 3" ψ 4" ψ 5 ψ 4" ψ 6" 4 ψ6 ψ 1" ψ 6 ψ 1" ψ 2" ψ 6 ψ 2" ψ 3" ψ 6 ψ 3" ψ 4" ψ 6 ψ 4" ψ 5" ψ 6 ψ 5" ψ 7" 4 ψ7 ψ 1" ψ 7 ψ 1" ψ 2" ψ 7 ψ 2" ψ 3" ψ 7 ψ 3" ψ 4" ψ 7 ψ 4" ψ 5" ψ 7 ψ 5" ψ 6" ψ 7 ψ 6" ψ 8" 4 ψ8 ψ 1" ψ 8 ψ 1" ψ 2" ψ 8 ψ 2" ψ 3" ψ 8 ψ 3" ψ 4" ψ 8 ψ 4" ψ 5" ψ 8 ψ 5" ψ 6" ψ 8 ψ 6" ψ 7" ψ 8 ψ 7" ψ 9" 4 ψ9 ψ 1" ψ 9 ψ 1" ψ 2" ψ 9 ψ 2" ψ 3" ψ 9 ψ 3" ψ 4" ψ 9 ψ 4" ψ 5" ψ 9 ψ 5" ψ 6" ψ 9 ψ 6" ψ 7" ψ 9 ψ 7" オリジナルは行列ベクトル積 ψ 8" ψ 9 ψ 8" 2016年5月12日 計算科学技術 特論B M EOT UP ベクトル積を行列積に変換 (DGEMV) 再帰分割法 (DGEMM) • • 依存関係のある三角部とない四角部にブロック化して計算 再帰的にブロック化することで四角部を多く確保 2016年5月12日 計算科学技術 特論B ※SDも同様に行列積化が可能 18 2016年5月12日 計算科学技術 特論B 3 8F7 B e gB A+ tm 63 l MWW P_O p m x 6 UYU Ud yu O(N2) ψ m H KS ψ n . lk ψn ψn M EOT UP B b 6+ gB A lk 8 9 8 9 H 8 9 8 9 8 9 H 8 9 l H m ,n = ψ m H KS ψ n O(N3) s B b 6+ gB A, lk m 8=5 z v M9 n −1 ' n ψ = ψ n − ∑ψ m ψ m ψ n m =1 P c P & $ $ $ % H N ×N N DZ H #& # & # !$ ! ! $! ! c = ε !$ n ! $ cn ! !$ ! $ ! "% " % " ψ n' ( r ) = ∑ cn ,mψ m ( r ) m =1 lk O(N3) B 6, gB A, O(N3) B b 6+ gB A, O(N3) i C Pc 8 9 P %8 9 H EOMWM MOV lk lk 8 9 8 9 H 8 9 D=77 1%))) 1b1b1 * E7: l l *+)b*+)b*+) × p E7: *)) )- × u 6OM 22 / B A, 8=5 ,) . B A, 8 9 8F7 +0 - B A+ ,1 / B A, 8 9 E7: 8B E7: */%))) = YP(= O *)) M EOT UP UbUYS% 3 D P_O % = YP(= O CE= C8E 9H8 6OM 3 6OM D P_O = = YP(= O YP(= O MWW P_O ( CE= (67E9F 5WW P_O ,* E7: z m 9A88B 2016年5月12日 計算科学技術 特論B F+ E 2) 6 5E 8 9 %8 9 H F _V_NM C = M UOT+ PU_ 3*+1% +./% .*+% *)+-% +)-1 Si4096(格子:96x96x96,バンド:8192) 実行時間,T2K-Tsukuba Si4,096(格子:96x96x96,バンド:8192) 速度向上率,T2K-Tsukuba 18.0 1000.0 SCF DIAG GS DTCG 900.0 Ideal SCF DIAG GS DTCG 16.0 14.0 12.0 700.0 速度向上率 実行時間(秒) 800.0 600.0 500.0 400.0 10.0 8.0 6.0 300.0 200.0 4.0 100.0 2.0 0.0 0 500 *+1 .*+ -b-b1 1b1b1 21 m 1,000 1,500 並列数 *)+1b1b*/ 2016年5月12日 計算科学技術 特論B 2,000 2,500 0.0 0 +)-1 1b*/b*/ 500 *+1 .*+ -b-b1 1b1b1 1,000 1,500 並列数 *)+1b1b*/ 2,000 +)-1 1b*/b*/ 2,500 Si4 0 9 6(格子:9 6x9 6x9 6,バンド:8 1 9 2) 演算時間と通信時間,T 2 K -Tsukuba 1000.0 通信(S :大域) 通信(S :隣接) 演算 900.0 SCF 800.0 実行時間(秒) 700.0 DIAG 600.0 500.0 400.0 GS 300.0 200.0 100.0 DTCG 128 256 512 1,024 2,048 128 256 512 1,024 2,048 128 256 512 1,024 2,048 128 256 512 1,024 2,048 0.0 SCF DIAG 並列数 GS DTCG fC8E 9H8 m 2016年5月12日 計算科学技術 特論B 1 2 3 4 5 6 2016年5月12日 計算科学技術 特論B 24 並列数 128 256 512 1,024 2,048 128 256 512 1,024 2,048 128 256 512 1,024 2,048 128 256 512 1,024 2,048 演算 通信(S :隣接) 通信(S :大域) 749.472 77.435 60.057 363.281 52.715 59.732 172.218 38.797 64.133 87.206 33.345 73.494 59.728 28.637 108.404 322.212 10.005 14.931 181.852 6.831 9.690 86.939 4.492 16.882 46.392 4.061 17.659 31.597 4.229 34.866 148.354 0.000 9.629 73.964 0.000 10.070 37.479 0.000 9.969 16.900 0.000 8.008 11.451 0.000 10.687 278.906 67.429 35.497 107.465 45.884 39.972 47.800 34.305 37.283 23.914 29.284 47.827 16.679 24.408 62.851 3 8F7 B e gB A+ tm p 6 UYU Ud l MWW P_O 通信時間増大 演算時間と逆転 . m 行列ベクトル積 x yu O(N2) ψ m H KS ψ n 63 性能は悪い 並列度の不足 lk ψn ψn M EOT UP B b 6+ gB A H m ,n = ψ m H KS ψ n O(N3) lk 通信時間減少せず l 8 9 s m 演算時間と同程度 8 9 H 8 9 行列積化で良好 並列度の不足 8=5 z v n −1 M9 ψ n' = ψ n − ∑ψ m ψ m ψ n m =1 P c P & $ $ $ % H N ×N #& # & # !$ ! ! $! ! !$ c n ! = ε $ c n ! !$ ! $ ! "% " % " N DZ H ψ n' ( r ) = ∑ cn ,mψ m ( r ) m =1 RSDFT • 実空間差分法 • ベクトルの内積計算 が基本 • 空間並列 B b 6+ gB A, O(N3) B 6, gB A, O(N3) B b 6+ gB A, O(N3) lk 通信時間増大 lk 8 9 行列積化で良好 演算時間と同程度 並列度の不足 i C Pc 8 9 P %8 9 H Scalapackの性能 EOMWM MOV Scalapackのスケー ラビリティが悪い lk 8 9 H 8 9 が悪い 8 9 lk 8 9 H 8 9 行列積化で良好 計算コアの最適化 • 行列積化 ターゲット計算機:PACS-CS,T2K-Tsukuba スレッド並列の実装 ターゲット計算機:PACS-CS,T2K-Tsukuba 超並列向けの実装 • バンド並列の拡張 EIGENライブラリ※の適用 ※高速固有値ライブラリ • Imamuraelal.SNA+MC2010(2010) ターゲット計算機:Kcomputer 2016年5月12日 計算科学技術 特論B 26 Hϕ i (r) = ε iϕ i (r) φ : 電子軌道(=波動関数) i i:電子準位(=エネルギーバンド) r:空間離散点(=空間格子) iはエネルギーバンド量子数 ML2 iについての依存関係はない 空間(S)に加えエネルギーバンド(B) の並列を実装 万を超える並列度を確保 3 2016年5月12日 計算科学技術 特論B 27 • • 10 • • • × • • 2016年5月12日 計算科学技術 特論B 28 - Gram-Schmidt 0 1 2 0 ? 1 ? 2 4 3 → 5 3 7 6 - 8 3 9 6 10 (1) 三角部の計算 (2) 計算した値を四角部に転送(バンド方向の各プロセッサに分配) (3) 四角部を並列に計算 HPCS2012 2016年5月12日 計算科学技術 特論B 29 - 通信の見積りと効果の予測 E 6 h h t ✓5 D98G79 ➢ M EOT UP 3 ➢8F7 3 ✓D98G79 3 ➢8=5 M9 ✓675EF3 ➢8=5 DZ H ➢ M EOT UP 3 5F 9DH ✓5 ➢8=5 ➢ M EOT UP ✓ ✓ ✓ 2016年5月12日 計算科学技術 特論B 367E9F CE= ux 8=5 M9 x - 通信の見積りと効果の予測 M EOT UP 3 ULMWWSM T UL MW1 6( ULMWW P_O UL MW1 A6 A6 ULMWW P_O UL MW1 A6 ULMWW P_O UL MW1 * ULNOM UL MW1 MWWSM T = * A6 * * * 6(A6 =Y WZS A6 A6 6(A6 =Y WZS A6 A6 * ) A6 6(A6 ( 6(A6 ( ( 6(A6 ( (A6 * 6(A6 ( ( 6 6 UL MW1 6 6 UOLNOM UL MW1 6E=K9 A6E=K9 ULU UL MW1 W MLY YP U MYV 6 / t 6( 6 ( UL MW1 W MLY YP U MYV 6 / t 6( 6 ( YP(U O P c * 6( 6 ULU YP 6( 6E=K9 6( 6 YP ( P 6 / 6( 6 ( UL MW1 P 6 / 6( 6 ( 6( 6 6 6(A6E=K9 ( UL MW1 ULaMU MWW v ( l ULU O 3 6( 6 * ULaMU MWW 63 A6 6E=K93 6b 6 6(A6 6(A6 UL MW1 ULU O 67E9F * 6(A6 ( (A6 * 6( EOMWM MOV P CE= A6 * * UL MW1 UL P_O 8=5 * A6 *3 3 UY 6E=K9%A6E=K9 v P3 ( 6E=K93 6b 6 W MLY YP3 - 通信の見積りと効果の予測 - 8F7 l l OZYPLOS 67E9F ULMWW P_O UL MW1 6LP 6( 6LP( ULMWW P_O UL MW1 6LP 6( 6LP( ULMWW P_O UL MW1 6LP 6( 6LP OS( ULMWW P_O UL MW1 6LP / 6( 6LP OS( ULMWW P_O UL MW1 6LP 6( 6LP OS( ULMWW P_O UL MW1 6 ULMWW P_O UL MW1 6LP 6( 6LP OS( ULU UL MW1 P 6LP / 6( 6LP OS( UL MW1 P 6LP / 6( 6LP OS( YP ULU O + ULaMU MWW CE= ULU 6( 6LP YP UL MW1 W MLY 6LP YP U MYV / UL MW1 W MLY 6LP YP U MYV / l 63 l ULU O l l ULaMU MWW l 67E9F ULU ULU O l ULaMU MWW 6LP 2016年5月12日 計算科学技術 特論B OS * ( OS( t 6( 6LP t 6( 6LP 6( 6LP YP l OS * ( OS * ( UL MW1 P 6LP / 6( 6LP OS * ( UL MW1 P 6LP / 6( 6LP OS * ( 6( 6LP P3 , W MLY YP3 OS * ( - 効果の確認 - I MV EOMWUYS ( v *+b*+b*+ 1 u tw ll F+ +%-)) * .*+ -1b-1b-1 *2%+)) .*+ -b-b-b1 + *%))) /)b/)b/) *2%+)) *%))) .b.b.b1 , *%0+1 0+b0+b0+ *2%+)) *%0+1 /b/b/b1 - -%)2/ 2/b2/b2/ *2%+)) -%)2/ 1b1b1b1 . 1%))) *+)b*+)b*+) *2%+)) 1%))) *)b*)b*)b1 F _V_NM 33 2016年5月12日 計算科学技術 特論B - 効果の確認 - I MV EOMWUYS S i4 0 9 6(格子:9 6x9 6x9 6,バンド:8 1 9 2) GS, Weak Scaling, T2K-Tsukuba Si4096(格子:96x96x96,バンド:8192) DTCG, Weak Scaling, T2K-Tsukuba 45.0 実行時間(秒) 50.0 40.0 30.0 20.0 DTCG 演算 通信(S:隣接) 通信(S:大域) 通信(B:大域) 40.0 35.0 実行時間(秒) GS 演算 通信(S :大域) 通信(B :大域) 60.0 30.0 25.0 20.0 15.0 10.0 10.0 5.0 0.0 0 1,000 2,000 3,000 並列数 4,000 0.0 5,000 0 20.0 15.0 10.0 5.0 4,000 5,000 RotV 演算 通信(S:大域) 30.0 実行時間(秒) 実行時間(秒) 25.0 3,000 並列数 35.0 MatE 演算 通信(S:隣接) 通信(S:大域) 30.0 2,000 Si4096(格子:96x96x96,バンド:8192) RoTV/DIAG, Weak Scaling, T2K-Tsukuba Si4096(格子:96x96x96,バンド:8192) MatE/DIAG, Weak Scaling, T2K-Tsukuba 35.0 1,000 25.0 20.0 15.0 10.0 5.0 0.0 0.0 0 1,000 2,000 3,000 並列数 4,000 5,000 0 1,000 2,000 3,000 並列数 4,000 5,000 -Tofu 空間並列 - 空間並列+バンド並列 バンド(1:On) Tofuネットワークへのマッピング 各バンドグループをサブメッシュ/ トーラス・ネットワークにマッピング バンド(1:Omax) バンド(On+1:Om) orbital z Z space Y y バンド(Om+1:Omax) x CPUspace X • • 2016年5月12日 計算科学技術 特論B マッピング・ルール 通信の最適化 サブメッシュ/トーラス内で通信が閉じられる HPCS2012 35 -Gram-Schmidt lgl llllllllll C= 120.0 実行時間(秒) 90.0 MPI_Bcast FZR_ s t u • 原子数: 19,848 • 格子数: 320x320x120 • 軌道数: 41,472 • トータルプロセス数: 12,288 mll ✓ 空間並列: 2,048(32x32x2) ✓ バンド並列: 6 60.0 MPI_Allreduce • 32x32x12のトーラスにマッピング 30.0 0.0 マッピングなし 2016年5月12日 計算科学技術 特論B 最適マッピング 36 -二軸並列の - SiNW,19,848原子,格子数:320x320x120,バンド数:41,472 トータル並列プロセス数は12,288で固定 120.0 Space +Orbital Space Space +Orbital Space Space +Orbital Space Space +Orbital Space TimeperSCF(sec.) in grid points are kept within a part and -79% do not affect the-78%Data-1 onWait/orbital 12288, 24576, 26864 and 73728 cores (cases 1, 2, 3 Globalcommunication/orbital communications within the other part. The number of parts is and 4, respectively). The number of parallel tasks is the product of 100.0 Globalcommunication/space taken up to the number of parallelisms in orbitals. We used this the number of parallel tasks in space and that in orbitals. The Adjacentcommunication/space mapping in the experiments. number of parallel tasks in space was fixed to 1,536, and the Computation 80.0 number of parallel tasks in orbitals was varied from 1, 2, 3 and 6. Therefore, the number of parallel tasks amounted to 1536, 3072, 5. PERFORMANCE ANALYSIS -78% 4608, 60.0 in grid points are kept within a part and do not affect or the9216. Data-1 on 12288, 24576, 26864 and 73728 cores (cases 1, 2, 3 The scalability and efficiency of the modified code which is communications within the part. above The number of parts is and 4, respectively). The number of parallel tasks is the product of parallelized in terms of grid points and orbitals as other described Figure 6 plots the computation time and communication time for taken up to the number of parallelisms in orbitals. We used this the number of parallel tasks in space and that in orbitals. The were evaluated with mapping a simulation model named “Data-1,” 40.0 GS, CG, MatE in SD of (MatE/SD) andinRotV SDfixed (RotV/SD) in the experiments. number parallel tasks space in was to 1,536, . and the -79% depicting a SiNW with 19,848 atoms. The number of grid points The horizontal axis in each figure stands for the number of cores. number of parallel tasks in orbitals was varied from 1, 2, 3 and 6. and orbitals were 320x320x120 and 40,992, respectively. We The theoreticalTherefore, computation time of forparallel cases tasks 2, 3amounted and 4 isto also the number 1536, 3072, 5. PERFORMANCE ANALYSIS 20.0 used an MPI library tuned the Tofuand network to obtain 4608, or 9216. time of case1 is divided by 2, 3 and the computation The to scalability efficiency of the higher modified codeploted whichwhere is performance. We also made an appropriate taskpoints mapping to theas described 6, respectively. They show the and parallelization fortime for parallelized in terms of grid and orbitals above Figure 6 plotsthe thescalability computationintime communication 空間分割:12,288 Tofu network. 0.0 were evaluated with a simulation model named the “Data-1,” orbitals. The timeand becames to GS,measured CG, MatE computation in SD (MatE/SD) RotV incloser SD (RotV/SD) . depicting a SiNW with 19,848 atoms. The number of grid points The horizontaltime axis in figure stands for the number the theoretical computation as each the number of parallel tasks inof cores. 5.1 Performanceand of orbitals DGEMM were 320x320x120 and 40,992, respectively. The theoretical time and for cases 2, 3 and orbitalsWeincreased, and GS, computation CG, 空間分割:2,048 MatE/SD RotV/SD were4 is also an MPI library tuned to the network to obtain higher The principal kernel of used the modified RSDFT code is aTofu DGEMM ploted where the computation time of case1 is divided by 2, 3 and found to be well-parallelized in orbitals. also made an appropriate 6, respectively. They show the scalability in the parallelization for routine of the BLAS. Weperformance. optimized We the DGEMM routine so as task to mapping to the Tofu such network. The communication time with to computation the parallel time tasksbecames in spacecloser to the orbitals. The respect measured バンド分割:6 utilize the special features as SIMD instructions and the the theoretical computation time as the number of parallel tasks in decreased in proportion to the number of parallel tasks in orbitals, sector cache and hardware synchronization capability of 5.1 Performance of DGEMM GSbarrier MatE/SD RotV/SD CG orbitals of increased, and GS,functions CG, MatE/SD and RotV/SD because the number calls of MPI for global and were the SPARC64TM VIIIfx The effectively. The sustained performance of principal kernel 行列生成 of the modified RSDFT 回転 code is a DGEMM found to be well-parallelized ininorbitals. adjacent communications decreased proportion to the the BLAS. We optimized the DGEMM DGEMM on a compute routine node isof123.7 giga-flops, or 96.6% of the routine so as to The communication timetasks with respect to theOn parallel tasks in space reciprocal of the number of parallel in orbitals. the other utilize the special features such as SIMD instructions and the peak performance. In particular, we found that the computation HPCS2012 decreased in proportiontime to thefor number paralleltasks tasks in sector cache and hardware capability of global hand, the communication the of parallel inorbitals, time as a result of keeping the block data onbarrier the synchronization L1 cache TM 2016年5月12日 計算科学技術 特論B because the number of calls of MPI functions for global and the SPARC64 VIIIfx effectively. The sustained performance of orbitals was supposed to increase as the number of parallel tasks manually decreased by 12% compared with the computation time adjacent communications decreased in proportion to the DGEMM on a compute node is 123.7 giga-flops, or 96.6% of the increased. in orbitals number MPI processes requiring for the usual data replacement operationsInofparticular, the L1 cache. reciprocalThe of the number of of parallel tasks in orbitals. On the other peak performance. we foundThis that the computation communicationshand, for the the global parallel tasks in orbitals, however, wastasks in communication time for the parallel DGEMM tuned for the K computer was also used for the time as a result of keeping the block data on the L1 cache actually to a was relatively small numberas of nodes, orbitals supposed to increase thecompute number of parallel tasks LINPACK benchmark program. manually decreased by 12% compared with the computation timerestricted in orbitals increased. numbercommunications of MPI processesofrequiring and therefore, the wall clock time The for global for the usual data replacement operations of the L1 cache. This communications for small. the parallel in we orbitals, however, was 5.2 Scalability DGEMM tuned for the K computer was also usedthefor parallel in orbitals was This tasks means succeeded the tasks actually to a relatively small of compute nodes, LINPACK benchmark We measured the computation time for theprogram. SCF iterations with in decreasing time for restricted global communication by number the combination 大域通信時間を大幅に削減 5.2 Scalability (a) We measured the computation time for the SCF iterations with 400.0 160.0 theoretical computation (a) 160.0 global/orbital theoretical computation 120.0 wait/orbital computation global/space Time per CG (sec.) 300.0 100.0 wait/orbital 200.0 80.0 40.0 80.0 40.0 100.0 0.0 0.0 (c) 0.040,000 Number0of cores 60,000 20,000 40,000 Number of cores (c) 200.0 Time per MatE/SD (sec.) 100.0 60,000 theoretical computation 150.0 (d) 40000 80000 60000 300.0 theoretical computation computation computation computation adjacent/space adjacent/space adjacent/space adjacent/space global/space global/space global/space global/space global/orbital global/orbital global/orbital global/orbital 50.0 200.0 100.0 80000 Number of cores theoretical computation 100.0 theoretical computation 200.0 100.0 0.0 0 20,000 60000 20000 Number of cores (d) 0.0 0 40000 0 computation 50.0 0.0 20000 0.0 80,000 300.0 200.0 150.0 0 80,000 Time per RotV/SD(sec.) 20,000 Time per RotV/SD(sec.) 0 Time per MatE/SD (sec.) adjacent/space global/space global/orbital 120.0 global/orbital Time per CG (sec.) Time per GS (sec.) 300.0 200.0 theoretical computation computation adjacent/space theoretical computation global/space computation global/orbital (b) computation global/space 400.0 Time per GS (sec.) 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 (b) 20,000 40,000 Number of cores 40,000 Number of cores 60,000 80,000 60,000 0 80,000 0.0 0 20,000 20,000 40,000 Number of cores 40,000 Number of cores 60,000 60,000 80,000 Figure 6. Computation and communication time of (a) GS, (b) CG, (c) MatE/SD and (d) RotV/SD for different numbers of cores. Figure 6.特論B Computation and communication time of (a) GS, (b) CG, (c) MatE/SD and (d) 2016年5月12日 計算科学技術 RotV/SD for different numbers of cores. 80,000 Hasegawa et al. 13 Table 2. Distribution of computational costs for an iteration of the SCF calculation of the modified code. Procedure block SCF SD MatE/SD EigenSolve/SD RotV/SD CG GS Execution time (s) Computation time (s) 2903.10 1796.97 525.33 492.56 779.08 159.97 946.16 1993.89 1281.44 363.18 240.66 677.60 43.28 669.17 Communication time (s) Adjacent/grids Global/grids Global/orbitals 61.73 823.02 12.57 13.90 497.36 4.27 13.90 143.98 4.27 – 251.90 – – 101.48 – 47.83 68.85 0.01 http://hpc.sagepub.com/ – 256.81 8.29 Performance Wait/orbitals (PFLOPS/%) 11.89 – – – – – 11.89 5.48/51.67 5.32/50.17 6.15/57.93 0.01/1.03 8.14/76.70 0.06/0.60 6.70/63.10 International Journal of High Performance Computing Applications The test model was a SiNW with 107,292 atoms. The numbers of grids and orbitals were 576 ! 576 ! 180, and 230,400, respectively. The numbers of parallel tasks in grids and orbitals were 27,648 and three, respectively, using 82,944 compute nodes. Each parallel task had 2160 grids and 76,800 orbitals. Article The International Journal of High Performance Computing Applications 1–21 The Author(s) 2013 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/1094342013508163 hpc.sagepub.com Performance Performance evaluation of ultra-largeevaluation of ultra-largescale first-principles electronic structure on the ª confinement becomes prominent. Thecalculation quantum code effects, K computer scale first-principles electronic structure depend on theAtsushi crystallographic of theTaisuke nano- Boku, Yukihiro Jun-Ichi Iwata, Miwako Tsuji, which Daisuke Takahashi, Oshiyama, directions Kazuo Minami, calculation code on Hasegawa, the K computer Hikaru Inoue, Yoshito Kitazawa, Ikuo Mitsuo Yokokawa wire axes andMiyoshi on the and cross-sectional shapes of the nanowires, Yukihiro Hasegawa et al., International Journal of High Performance Computing Applications published online 17 October 2013 result in substantial modifications to the energy-band DOI: 10.1177/1094342013508163 Yukihiro Hasegawa1, Jun-Ichi Iwata2, Miwako Tsuji1, structures and the transport characteristics of SiNW FETs. 3 2 1 Daisuke Takahashi , Atsushi Oshiyama , Kazuo Minami , The online version of this article can be found at: However, knowledge of the effect of the structural mor, Taisuke Boku3, Hikaru Inoue4, Yoshito Kitazawa5http://hpc.sagepub.com/content/early/2013/10/16/1094342013508163 Ikuo Miyoshi6 and Mitsuo Yokokawa7,1 phology on the energy bands of SiNWs is lacking. In addition, actual nanowires have side-wall roughness. The 2016年5月12日 計算科学技術 特論B Published by: imperfections on the energy bands are effects of such Abstract Silicon nanowires are potentially useful in next-generation field-effect transistors, and it is important to clarify the electron unknown. The advent of reliable first-principles calculastates of silicon nanowires to know the behavior of new devices. Computer simulations are promising tools for calculating http://www.sagepublications.com electron states. Real-space density functional theory (RSDFT) code performs first-principles electronic structure calculations would provide a firm theoretical framework for the tions. To obtain higher performance, we applied various optimization techniques to the code: multi-level parallelization, load balance management, sub-mesh/torus allocation, and a message-passing interface library tuned for the design K computer. of suitable SiNWs of 10,000–100,000 atoms for We measured and evaluated the performance of the modified RSDFT code on the K computer. A 5.48 petaflops (PFLOPS) Additional services and information for International Journal of High Performance Computing be imporfound at: Our RSDFT code enables usApplications to achievecan this sustained performance was measured for an iteration of a self-consistent field calculation for a 107,292-atomFETs. Si nanowire simulation using 82,944 compute nodes, which is 51.67% of the K computer’s peak performance of 10.62 PFLOPS. This tant task using the K computer. scale of simulation enables analysis of the behavior of a silicon nanowire with a diameter of Email 10–20 nm. Alerts: http://hpc.sagepub.com/cgi/alerts We have performed extensive electronic-structure calKeywords Subscriptions: http://hpc.sagepub.com/subscriptions culations for [110]- and [100]-SiNWs with various crossK computer, Tofu network, next-generation supercomputer, real-space density functional theory, RSDFT, self-consistent electron states, silicon nanowire, petaflops, PFLOPS, collective communication Reprints: http://www.sagepub.com/journalsReprints.nav sectional shapes using RSDFT with local density approxiFigure 13. Plot of volume versus total energy for crystalline mation (LDA) (Perdew and Zunger, 1981). Figure 14 supercomputer (http://www.open-supercomputer.org/; NaPermissions: http://www.sagepub.com/journalsPermissions.nav 1. Introduction silicon of 10,648 (squares) and 21,952 (triangles) kashima, 2008),atoms. by incorporating matrix–matrix multipli-examples of cross-sectional views of the obtained shows Computer simulations are essential when clarifying and precations for Gram-Schmidt orthogonalization with a high dicting the properties of materials that have prospects for geometry-optimized [110]-SiNWs. The cross-sectional cache-hit ratio and thread parallelization with OpenMP practical applications. In particular, first-principles electronic directives for multicore processors (Iwata et al., 2010; Tsuji 3 3 structure-times calculationsand based on functional theory (11 14density -times supercells correspond to the dimension is approximately 3–6 nm. The sidewalls of the >> OnlineFirst Version of Record - Oct 17, 2013 (DFT) have been performed on a variety of materials using 10,648and 21,952-atom models, respectively). The comwires are terminated by hydrogen atoms, reflecting certain diverse software implementations on parallel computers. RIKEN Advanced Institute for Computational Science, Kobe, Japan What is This? Real-space density functional theory (RSDFT) code Department of Applied Physics, School of Engineering, The University of experimental situations on one hand and focusing on essenputed equilibrium lattice constants for the two supercell (Iwata et al., 2010) is a simulation technique used to per- Tokyo, Tokyo, Japan form first-principles electronic structure calculations. ‘Real Center for Computational Sciences, University of Tsukuba, Tsukuba, tial features of the energy bands on the other. Figures 14(b) models correspond to the lattice constant, 5.39 Å, of the space’ means that three-dimensional physical coordinates Japan Computational Science and Engineering Solution Division, Technical are discretized, and the wave functions, electronagrees density, with and (d) show ellipses that extend along different crystalloprimitive unit cell, which the experimental Computing Solution Unit, Fujitsu Ltd, Chiba, Japan and potential field are calculated at the resulting discrete IT Solution Unit, CAE Simulation Department, Fujitsu Systems East Ltd, result (5.43 withinof 1%. graphic directions. We call the ellipses that have their long lattice points or grids.Å) One advantage this method is that Nagano, Japan it is suitable for parallel computations. In fact, the HamiltoPA Project, Next Generation Technical Computing Unit, Fujitsu Ltd, Prediction offormulation electronic axes along the [001] and the [1-10] directions the [001]nian matrix of the real-space is sparse, andstructures, fast Kawasaki, Japan or correlation Graduate School of System Informatics, Kobe University, Kobe, Japan Fourier transform (FFT), which usually requires global combetween the electronic structure and a specific atomic ellipses and the [1-10]-ellipses, respectively. Similarly, the munication traversing all compute nodes of a parallel com- Received June 13, 2013; Accepted for publication September 9, 2013. puter, is unnecessary Hamiltonian matrix operations. application of DFT. To structure, isforanother important dumbbells shown in Figure 14(c) and (e) are called the Corresponding author: RSDFT code has been parallelized to run on parallel Yukihiro Hasegawa, RIKEN Advanced Institute for Computational Science, investigate possibility, we applied the present RSDFT computers, like thethis PACS-CS (http://www.rccp.tsukuba. [001]-dumbbells and [1-10]-dumbbells, respectively. Fig7-1-6, Minatojima-minami-machi, Chuo-ku, Kobe, Hyogo 650-0047, Japan. ac.jp/PACS-CS/; Boku et al., 2006) and the T2K open [email protected] scheme to SiNWs. Si nanowire Email: field-effect transistors ure 15 shows the calculated energy bands for each SiNW (SiNW FETs) are expected to be boosters in post-scaling shown in Figure 14. In all cases, the conduction band minsemiconductor technology (http://www.itrs.net/). Clear imum is located at the zone center !. The second minimum, scaling of short-channel effects vs. nanowires with fixed which is higher than the !-point minimum by "Eshift, is gate lengths has been observed (Bangsaruntip et al., located at some point between ! and the zone boundary. 2009). Reducing the dimensions of SiNWs has been The "Eshift is a deciding factor in determining the number shown to improve short-channel control. Furthermore, of conduction channels and is therefore important when SiNW FETs with channel dimensions of 5.0 ! 6.3 nm designing FETs. We have found that the "Eshift of the have been fabricated (Bangsaruntip et al., 2009). A recent [001]-ellipse and dumbbell are smaller than those of the experiment showed that the optimal dimension of the [1-10]-ellipse and dumbbell, and that the "Eshift of cross-section of a SiNW FET is approximately 10 nm the [001]-ellipse is smaller than that of the circle. These (Chen et al., 2010). At such a small scale, quantum results unequivocally show that, with respect to the number 2016年5月12日 計算科学技術 特論B 1 2 3 4 5 6 7 Downloaded from hpc.sagepub.com at RIKEN on April 30, 2014 Downloaded from hpc.sagepub.com at RIKEN on April 30, 2014 Downloaded from hpc.sagepub.com at RIKEN on April 30, 2014 ● ● ● 2016年5月12日 計算科学技術 特論B Kohn-Sham方程式 Hϕ i (r) = ε iϕ i (r) 波数:Gによる展開 Hϕ ik (G) = ε iϕ ik (G) φik: 電子軌道(=波動関数) i:電子準位(=エネルギーバンド量子数) G:波数格子 k:k点 2016年5月12日 計算科学技術 特論B 42 解くべき方程式 H|ψ> =εn|ψ> |ψ> i : Δi+1 <= F( H|ψ>i , |ψ> i ) |ψ>i+1 <= |ψ>i + Δ i+1 |ψ>i+1 H|ψ> |ψ> |ψ> |ψ> |ψ>i |ψ> <ψ|ψ>i+1 : H : 2016年5月12日 計算科学技術 特論B Hϕ ik (G) = ε iϕ ik (G) iはエネルギーバンド量子数 基本的にエネルギーバンドについて並列化されてい る 一部,波数:Gについて並列化されている G G並列の前にエネルギーバンド並列されている波動 関数をG並列可能なようにトランスバース転送が発生 G並列後にG並列されている波動関数をエネルギーバ i i 2016年5月12日 計算科学技術 特論B ンド並列に戻すためなトランスバース転送が発生 このトランスバース転送のコストが大 G ■ O(N3) O(N2logN) O(N3) 2016年5月12日 計算科学技術 特論B BLAS Level3 LAS! Level3 ! ■ ! BLAS Level3 BLAS Level3 ( ( ( 2) ( 2) 14 [s] 14 [s] 12 12 10 10 8 8 6 6 4 4 2 2 0 origina 23axis: 23axis: 0 16 46 original:0calc. 23axis:0calc. 23axis:0comm. 32 16 64 32 2016年5月12日 計算科学技術 特論B 12864 1 ■ 第 4 章 性能最適化技術の応用 演算処理 通信処理 トランスバース処理 (通信以外) 2016年5月12日 計算科学技術 特論B Fig. 4.24: PHASE 区間 5 のスケーラピリティ 確保できないということである。良い並列効果を得るには、オリジナルコードを使 限り、ノードあたりの原子数をもっと大きくとる必要がある。区間 8 でも、区間 2 のハードウェアと現状のアプリの並列度のミスマッチが発生していると言える。 ■ 41 41 41 41 41 41 5637 41 41 0. 2 Fig. 4.25: PHASE 区間 8 のスケーラピリティ 2016年5月12日 計算科学技術 特論B 第 4 章 性能最適化技術の応用 対角化 対角化の処理を実施している区間 9 を HfSiO2、384 原子アモルファス系を入力データと ■ して、16 から 128対角化カーネル(区間9) 並列の低並列度でストロングスケールで測定した。スケーラビリティの データを Table 4.8 に示す。固有値計算の部分は並列度を上げても全く実行時間の短縮が見 HfSiO2 384原子,アモルファス系で測定. られない。これは並列対応のライブラリである ScaLAPACK でなく、逐次対応の LAPACK 固有値計算の部分は並列度を上げても全く実行時間の短縮が見られない . を使用していることが原因である。つまりこの部分は、完全に非並列部が残存しているこ ととなっている。逐次対応のLAPACK を使用していることが原因. 完全に非並列部が残存している. Table 4.8: PHASE 区間 9 のスケーラピリティ プロセス数 16 32 64 128 区間 区間 9 全体 固有値計算 (lapack-dsyev) ! ψi′ = Uij ψj ′ ψi 転置通信 高速化率 [-](全体) 高速化率 [-](固有値 計算含まず) 22.233 21.228 20.881 20.633 20.349 20.189 20.260 20.223 1.584 0.824 0.410 0.210 0.300 0.215 0.211 0.200 1.00 1.05 1.06 1.08 1.00 1.81 3.03 4.60 2016年5月12日 計算科学技術 特論B 4.3.8 PHASE の高性能化のための書き換え 並列軸の拡張 4.3.7 節の行列行列積に書き換え可能な処理と FFT を含む処理の項で、PHASE の高並 列上の問題点は、ハードウェアとアプリケーションの並列度のミスマッチであることを述 φik: 電子軌道(=波動関数) べた。PHASE は (4.4) の固有値方程式を解いている。オリジナルの PHASE では基本的 ik i ik i:電子準位(=エネルギーバンド量子数) に、固有値方程式の中のエネルギーバンド量子数を表す変数、 i について並列化が実施さ G:波数格子 れている。また一部分、例えば GramShumit 対角化処理については、波数 (G) による並 列化が実施されている。本来 PHASE の方程式は、すべての部分が完全に、エネルギーバ k:k点 ンド:i と波数:G について並列化可能である。 3.6 節に述べたように並列軸の拡張により、 ハードウェアとアプリケーションの並列度のミスマッチへの対処が可能になる。今回は並 エネルギーバンド(B)に加え波数(G) 列軸の拡張として、エネルギーバンドと波数について、すべての部分について 2 軸のフル 並列化を実施した。この並列軸の拡張により、 の並列を実装 1000 止まりであった並列度を数万レベルま で拡張し、原子数よりも大きな並列度で並列化することを目指した。 完全な2軸並列とする また 4.3.7 節で PHASE の高並列上の問題点は、非並列部の残存であることを示した。 RSDFTと同様の行列積化も実装 Fig. 4.24 に示した、原子数ループの先頭とエネルギーバンド並列部に挟まれた非並列部 は、波数に関する計算を実施している部分である。したがってエネルギーバンドに加えて、 波数 (G) による並列化を実施することにより、この残存している非並列部が並列化可能と なる。 Hϕ (G) = ε ϕ (G) 2016年5月12日 計算科学技術 特論B 50 ■二軸並列化 ■ 非並列部が波数で並列 化できる 2016年5月12日 計算科学技術 特論B ■二軸並列化 ■ • • • 2016年5月12日 計算科学技術 特論B 51 ■二軸並列化 ■ • バンド 波数 2016年5月12日 計算科学技術 特論B ■二軸並列化 ■ • • • • • 波数 バンド 2016年5月12日 計算科学技術 特論B ■二軸並列化 • この分野では小規模問題を短時間で計算したいという科学 的要求が高い. • バンド計算(エネルギー準位など):1万原子の1回SCF収束 で良い∼100SCF程度. • 構造緩和(MD)や反応経路探索:外側に原子核の緩和に関す るループ構造∼100step程度. • 10,000原子を10PFシステム(80,000ノード),また1,000 原子を10,000ノードで計算する事を目指せる. • ただし二軸並列はメリットとデメリットがあるため実施前 に効果が期待できるか詳細な評価を実施した. 2016年5月12日 計算科学技術 特論B [sec] [s] 300 250 200 [%] [sec] [%] [sec] [%] [sec] [%] [%] [%] [sec] [%] 5 512.7 100.0 2 23 5 512.7 100.0 2 23 [sec] [%] [sec] [%] 512.7 100.0 [%] 105.3 m_ES_F_transpose_r 105.3 20.5 0 m_ES_F_transpose_r 20.5 [%] 5 0 m_ES_F_transpose_r 105.3 20.5 5 512.7 100.0 2 23 m_ES_W_transpose_r 15.7 3.1 0 m_ES_W_transpose_r 15.7 3.1 m_ES_W_transpose_r 0 15.7 3.1 m_ES_F_transpose_r 105.3 20.5 0 WSW_t 15.2 3.0 36 WSW_t 15.2 3.0 36 WSW_t 15.2 3.0 m_ES_W_transpose_r 15.7 3.1 0 0.2 normalize_bp_and_psi_t 0.9 0.2 3 250.9 normalize_bp_and_psi_t 0.2 normalize_bp_and_psi_t 3 25 WSW_t 15.2 3.0 36 0.9 Gram/Schmidt3orthonormalization 49.2 9.6 W1SW2_t_r 49.2 9.6 5 46 normalize_bp_and_psi_t 0.9 0.2 3 25 49.2 W1SW2_t_r 9.6 W1SW2_t_r 5 46 modify_bp_and_psi_t_r 50.8 9.9 W1SW2_t_r modify_bp_and_psi_t_r 49.2 9.6 5 46 modify_bp_and_psi_t_r 50.8 9.9 4 45 50.8 9.9 4 45 W1SW2_t_r_block 162.0 31.6 modify_bp_and_psi_t_r 50.8 9.9 4 45 W1SW2_t_r_block 162.0 31.6 41 89 W1SW2_t_r_block 162.0 31.6 modify_bp_and_psi_t_r_block 41 89 other 96.2 18.8 W1SW2_t_r_block 162.0 31.6 41 89 M× M modify_bp_and_psi_t_r_block 96.2 18.8 74 65 modify_bp_and_psi_t_r_block 96.2 18.8 74 65 m_ES_W_transpose_back_r 14.1 2.8 modify_bp_and_psi_t_r_block 96.2 18.8 74 65 M× V 0.3 m_ES_W_transpose_back_r 14.1 2.8 01.3 m_ES_W_transpose_back_r 2.8 m_ES_F_transpose_back_r 0 m_ES_W_transpose_back_r 14.1 2.8 0 14.1 m_ES_F_transpose_back_r 1.3 0.3 0 m_ES_F_transpose_back_r 1.3 0.3 0 m_ES_F_transpose_back_r 1.3 0.3 0 150 行列積 100 行列ベクトル積 50 0 orinigal original+DGEMM 22axis+1threads Algorithm 22axis+8threads 2016年5月12日 計算科学技術 特論B [%] [%] 2 23 0 0 36 3 25 5 46 4 45 41 89 74 65 0 0 ■二軸並列化 ■ • 行列積化されたカーネルに(区間2)ついての結果. • HfSiO2 384原子アモルファス系のデータ. • 大幅な性能向上を達成. 通信部 20.0 20.0 18.0 18.0 16.0 14.0 12.0 ]ces[間時過経 ]ces[間時過経 16.0 14.0 12.0 10.0 8.0 6.0 4.0 2.0 0.0 16 32 64 14 [s] 主要演算ループ original:0calc. 23axis:0calc. 23axis:0comm. 12 10 8 10.0 8.0 6 6.0 4.0 4 2.0 0.0 2 128 0 プロセス数 16 プロセス数 2016年5月12日 計算科学技術 特論B 32 64 プロセス数 128 ■二軸並列化 ■FFT • FFTを含むカーネルに(区間8)ついての結果. • HfSiO2 384原子アモルファス系のデータ. • 性能向上を達成. 50.0 50.0 40.0 40.0 18 [s] ]ces[間時過経 ]ces[間時過経 14 30.0 30.0 12 20.0 10 10.0 10.0 8 0.0 0.0 20.0 16 32 64 original:0other original:0DGEMM original:0FFT 2:axis:0other 2:axis:0comm. 2:axis:0DGEMM 2:axis:0FFT 16 6 4 128 2 プロセス数 2016年5月12日 計算科学技術 特論B プロセス数 0 16 32 プロセス数 64 128 ■Scalapack分割数の固定 ■ • 対角化はエネルギーバンド数の元を持つ行列が対象 • 行列の大きさに比べて分割数が多すぎる • 分割数を16 16=256に固定 100 [s] 700[s] HfSiO4_3072 80 HfSiO4_6144 60 600 250[s] ScaLAPACK FFT BLAS 200 500 400 150 300 40 200 20 100 100 0 0 0 10 20 分割数 30 50 40 0 16x32 16x64 16x96 16x128 プロセス数(バンド×波数) 2016年5月12日 計算科学技術 特論B ■ ■ 16x160 ( Time%[s] 1400 DGEMM FFT ScaLAPACK 1200 1000 800 600 400 200 0 48 2016年5月12日 計算科学技術 特論B 96 192 384 768 1536 3072 6114 12288 Process 2016年5月12日 計算科学技術 特論B
© Copyright 2025 ExpyDoc