富士通HPCフォーラム.key

Transcription

1 & & RIKEN&ADVANCED&INSTITUTE&FOR&COMPUTATIONAL&SCIENCE はじめに LINPACK HPCGとは HPCGベンチマークプログラム HPCGのチューニングと性能 2

2 3 4

3 5 6 Jun Nov.2011 Jun Nov Jun Jun System country PFLOPS System country PFLOPS System country PFLOPS System country PFLOPS System country PFLOPS 1 K computer JPN 8.16 K computer JPN Sequoia USA Titan USA Tianhe-2 (CHN) Tianhe A CHN 2.57 Tianhe A CHN 2.57 K computer JPN Sequoia USA Titan USA Jaguar USA 1.76 Jaguar USA 1.76 Mira USA 8.16 K computer JPN Sequoia USA Nebulae CHN 1.27 Nebulae CHN 1.27 SuperMUC GER 2.90 Mira USA 8.16 K computer JPN TSUBAME2.0 JPN 1.19 TSUBAME2.0 JPN 1.19 Tianhe A CHN 2.57 JUQUEEN GER 4.14 Mira USA 8.59

4 7 High Performance LINPACK(HPL) TOP500& TOP500J.0Dongarra LINPACK & LINPACK TOP & LINPACK

5 9 N N( ) p(1 ) q(1 ) nb(1 ) nb 10 p q nb p q nb DTRSM () p q nb

6 11 q p nb q p nb DTRSM () DGEMM ( ) 成功 12 わかりやすさ システム最大性能との線形相関 問題点 (*)&HPCG WorkShop,25,March,2014,Bethesda, MD 実問題との乖離 発足して 20 年以上が経過 実際のアプリケーションで求められる性能要件との乖離 Stunt machine の開発を助長する 長すぎる実行時間 途中からの実行など実行時間短縮案を検討中

7 13 14 (*)&HPCG WorkShop,25,March,2014,Bethesda, MD ベンチマーク自体の紹介 HPCG(High Performance Conjugate Gradient) Version 2 の紹介 TOP500 との関係 置き換えるものではない TOP500 のリストに列を追加して HPCG スコアも表示するなどの案あり WS 参加者のなかでは将来的に置き換えるべきという意見もあった

8 15 左図のような 2 次元の領域で微分方程式を解くとする. 例として 5 点差分で差分化. 元の微分方程式は, 差分化により連立一次方程式を解く事に帰着. Ax = b f (x) = 1 2 ax2 bx f (x) = 1 (x,ax) (b,x) 2 f '(x) = ax b = 0 f (x) = 1 n a ik x i + 1 x k 2 2 i=1 n j=1 n f (x) = a ki x i b k = 0 x k i=1 a kj x j b k 16 f (x) = 1 (x,ax) (b,x) 2 CG 法アルゴリズムの基本形. この他にアルゴリズムの派生形あり. f(x) が最小値が得られるように xi の列を求める. ri は残差ベクトルであり互いに直交し一次独立である. 残差ベクトルは N 元の連立一次方程式には N 個しか存在しない. N 元の連立一次方程式は高々 N 回の反復で収束する. また A の固有値が縮重しているか密集していれば収束が速くなる性質を持つ.

9 17 M = U T U Ux = x U T b = b U T AU 1 = A A = U T AU 1 I Ax = b A に近い正値対称行列 :M を考えコレスキー分解する. このような x,b,a の置き換えを行う. 置き換えた A の性質は単位行列に近い. 置き換えた行列 A の固有値は 1 の周りに密集している 置き換えた行列とベクトルに対する連立一次方程式の収束は早いものと期待. CG 18 不完全コレスキー分解による前処理 CG 法アルゴリズムの基本形. 計算は行列ベクトル積 ベクトルスカラ積 ベクトルの和 内積 除算で構成される. 赤線で示した前処理部分は前進代入 後退代入で計算される. α k = (r i k (LL T ) 1 r i k )/(Ap i k p i k ) x i k+ 1 = x i k +α k p i k r i k+ 1 = r i k α k Ap i k β k = (r i k+ 1 (LL T ) 1 r i k+ 1 ) /(r i k (LL T ) 1 r i k ) p i k+ 1 = (LL T ) 1 r i k+1 + β k p i k

10 19 20 マルチグリッドも行列ベクトル積の計算が中心 f k = R k+1 k f k+1 f k 1 = R k k 1 f k f k 2 = R k 1 k 2 f k 1 f 1 = R 2 1 f 2 A k+1 u k+1 = r k+1 r k+1 = f k+1 A k+1 u k+1 u k+1 = R k k+1 u k A k u k = r k r k = f k A k u k u k = R k 1 ku k 1

11 A Processor Performance 21 High Performance Low Performance 22

12 Ax = b 連立一次方程式を前処理付 CG 法で解く際の性能をベ ンチマーク 前処理はVer.2.0 以降はMulti-Grid 法を採用 (Ver1.1まではGauss-Seidel 法 ) 係数行列 Aは有限要素法の行列 疎で対称 具体的には直方体領域を規則的に分割したもの 特徴 LINPACKと比して通信 メモリアクセスの比重が高い メモリアクセスはリストアクセスであり, より高いメモリ性能が必要 C++ で書かれているが便利なCとしてしか使っていない Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results 大まかに4つの部分に分けられる 前処理 メモリ確保などの初期化 問題設定, 行列の最適化 検証 ユーザーチューニングの正当性検証 雑多なデータ測定 測定 実行時間の予測 本番測定 (1 時間 ) 後処理 結果を整理し提出するデータを作成

13 25 Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results XYZ GenerateGeometry()& & Ax=b GenerateProblem()& & & SetupHalo() & Multi^Grid GenerateCoarseProblem()& & OptimizeProblem()& & & Ax=b ( ) P0 P )& 3 3 3& & 2 z y x y z & x & 2 & ELL / global_int_t (int) doulbe

14 MG 27 3 & 4 & 1, Generate, Geometry 2, Generate, Problem 3, SetupHalo 4, Generate, CoarseProblem XYZ & XYZ & etc & & etc & & & XYZ 1/2 1, Generate, Geometry 2, Generate, Problem 3, SetupHalo 4, Generate, CoarseProblem XYZ & XYZ & etc & & etc & & & 28 Problem&Setup&Phase & 2 ComputeCG_ref() ComputeCG() Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results カーネルのチューニングや行列の最適化が正しいかをチェック TestCG() チューニング版 CG 関数に, 強く対角優位な行列を与えて, 規程回数で収束するかを判定 ( 収束判定値 ) 前処理なしなら 12 回以内 前処理ありなら 2 回以内 TestSymmetry() OptimizeProblem() で最適化した行列が対 称になっているかチェック 条件は決め打ち

15 29 Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results 基準版の SPMV( 疎行列ベクトル積 ) と MG(Multi- Grid) 関数の実行時間を計測 ComputeSPMV_ref() 基準版の疎行列ベクトル積 ComputeMG_ref() 基準版の MG 10 この時の SPMV と MG の 1 回当りの平均時間を求める Report Result 部で OptimizeProblem() の時間と, この平均 時間との比を出力している 30 Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results 1 チューニング版 CGを1 回実行し, 要した時間を測定 CG() チューニング版 CG 反復上限は500 回 収束判定はCG_ref() を実行して求めた ref_tolerance(1.0e-3 程度 ) 1 回だけ実行 この時間を元に, 次のフェーズでCG() を何回実行するかを決定 ( ) HPCG

16 31 Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results 1 本番測定チューニング版 CG を複数回実行し測定 複数回 CG() チューニング版 CG 反復上限は先の1 回実行の反復数 ( 収束判定は 0 なので必ず指定回数反復 する ) 合計実行回数は, 先に計った1 回の実行時間を元に, 入力データで指定した時間となるよう決定 1 TestNorms() 再現性 安定性のチェック CGの実行毎にの値を記録 r 50 r しておき, ブレがないかチェック ^6 z = MG( A, r ) z = 0 A z = r Az = A z r coarse = Restrict( r Az) z coarse = MG( Acoarse, r z = z + Prolongate( z A z = r Relaxation&&&&&&&&&&&&&&&&&&&&&&&by&&Gauss^Seidel&&&&(Pre^smoothing) coarse coarse ) ) () & Relaxation&&&&&&&&&&&&&&&&&&&&&&&&by&&Gauss^Seidel&&&&(Post^smoothing) A z = r Relaxation&&&&&&&&&&&&&&&&&&&&&&&by&&Gauss^Seidel 32 ( ) & N N k+ 1 1 ( k+ ) ( )# 2 k i = % ri aij z j aij z j " aii $ j< i j> i Sweep&&&k =&1 N z 2 1 Sweep&&&k=N 1 ( ) & N N ( ) ( )# k+ 1 1 k+ k+ zi = % ri aij z j aij z j " aii $ j< i j> i MG& Restrict Z=0& Pre^smoothing& Residual& Restrict Prolongate )& MG& r A coarse z 1:2 Prolongate& MG& Post^Smoothing coarse Relax corase z

17 33 Problem&Setup&Phase Validation&Testing&Phase Reference&SpMV+MG&Timing&Phase Reference&CG&Timing&Phase Optimized&CG&Setup&Phase Optimized&CG&Timing&Phase Report&Results 提出する実行結果を出力 (HPCG-Benchmark-2.1- 年月日時刻.yaml) 並列数 スレッド数 全体問題サイズ プロセス形状 1ノード当りの問題サイズ マルチグリッド各階層での行列の行数と非ゼロ数 チューニングの正当性検証結果 再現性チェック結果 測定結果 ( 各カーネルと全体 ) Elapse 時間 理論的な浮動小数点演算数 GFLOPS 値 検証フェーズ で不合格があればFAILED, 正常ならVALID+GFLOPS 値が出力される 34

18 35 コードはAs Isそのまま ただし大並列への対応として一部の変数の型を変更 Geometry.hpp 中のtypedefされている型をintからlong longに変更 Make 方法 QUICKSTART.TXTに従って実行 テスト問題 デフォルトの =112 万節点 / プロセスを利用 測定時のCG 繰り返し回数は5 回 ( 本来は1 時間分で約 50 回 ) 36 反復上限は 50 回 CG 実行回数は 5 回に固定大並列で動かすよう一部変数の型を変更 15222GFLOPS&(0.36%) 9835GFLOPS&(0.23%) 8 並列 並列の 4096 倍で GFLOPS 値は 4054 倍 ( ベンチ結果 ) 4040 倍 ( プロファイラ ) 並列化効率 98% 8 / 8 並列から 並列まで並列数を換え測定 ベンチ結果 HPCGコードが自身で算出した ( 提出する ) GFLOPS 値 プロファイラプロファイラを利用して算出した GFLOPS 値 並列までの範囲ではウィークスケール性能は良好 3 ベンチ結果とプロファイラとで 1.5 倍程度の差がある 浮動小数点演算数カウント方式の違いが原因

19 37 プロファイラでの浮動小数点演算数 測定区間で発生する浮動小数点演算全てがカウント対象 除算等は複数の浮動小数点演算としてカウント HPCGでの浮動小数点演算数 理論的に必要な演算のうち, 主要なもののみをカウント しており, 無視しているものもある x = CG( A, b, x0, ε, max) x = SYMGS( A, b, x) SPMV() WAXPBY() DotProduct for(k=1; k<max && conv; k++){ MG() } DotProduct() WAXPBY() SPMV() V2.4 & for(i=0; i< ; i++){ for(j=0; j<nz[i]; j++){ col = A.mtxIndL[i][j] sum -= val[j] * x[col]; } sum += x[i] * Diag[i]; sum = sum / Diag[i]; } 38 4x4x4 8x8x8 6% 2x2x2 実行時のプロセス数指定は GenerateGeometry で決定した 3 次元形状に合わせた 3 プロファイラで計測した MPI のコストの, 対全コスト比率 通信のコストの割合は,512 並列以上ではサチってきており, 合計 6% 程度

20 39 ComputeDotProduct_ref&&46^79 ExchangeHalo_ref&&59^70 ComputeSPMV_ref&&59^70 ComputeSYMGS_ref&&54^103 プロファイラで採取した CG 内で呼ぶ手続きのコスト分布 上位 4 つで 98% 以上を占める (Allreduce, Isend, Irecv, Wait を含む ) チューニングはまず単体性能の向上 ComputeSYMGS と ComputeSPMV が主眼に 8192 () 40 for(int i=0; i<nrow; ++i){ matrixvalues[i] = new double[27]; mtxindl[i] = new int[27]; } new A.matrixValues[0] A.mtxIndL[nrow^1] A.mtxIndL[1] A.matrixValues[nrow^1] A.mtxIndL[0] A.matrixValues[1]

21 41 double* tmp1 = new double[nrow * 27]; double* tmp2 = new int[nrow * 27]; for(int i=0; i<nrow; ++i){ matrixvalues[i] = &(tmp1[i*27]); mtxindl[i] = &(tmp2[i*27]); } & A.matrixValues[0] A.matrixValues[1] A.matrixValues[nrow^1] A.mtxIndL[0] A.mtxIndL[1] A.mtxIndL[nrow^1] / 2.3 CG() から直接呼ばれる SPMV について, 実行時間と, キャッシュミス率, スループッ ト等を調べた MG() から呼ばれる分は除外 メモリ領域の連続化の効果によりキャッ シュミス率減少, スループット向上 実行時間は 1/2.3 に短縮

22 43 CG MGSPMV,SYMGS 16.2GFLOPS& 1.58% X GFLOPS& 3.82% 2.3GFLOPS& 0.22% X GFLOPS& 0.7% 2.6GFLOPS& 0.25% X GFLOPS& 0.82% HPCG 自身が測定した提出するデータから抜粋したもの SPMVの時間 CG() から呼ばれる分のみの時間 前処理 MGはSPMV,SYMGS 等々を含む時間 2 点の簡単なチューニングだけで3.2 倍の 高速化が得られた 行列のメモリ領域の連続化 SYMGS の BAK ループでのループ方向の反転 44 ComputeSYMGS_ref& &93^96 ComputeSYMGS_ref 74^77 ComputeSPMV_ref& &67^68 プロファイラで採取した各手続きのループ単位のコスト比率 ComputeSYMGS のループはマルチスレッド化されておらず逐次で動作 ここでは ComputeSYMGS のマルチスレッド化を実施

23 45 SYMGSは対称行列についてガウス ザイデル法を使用したマルチグリッド処理を行う. i 1 n ガウス ザイデル法 a ii x (m+1) (m+1) i = a ij x j a ij x (m) j + k j j=1 (m+1) 世代の解を得るために (m+1) 世代の下三角行列を使用して ( リカレンスの発生 ). j=1+1 (a) (c)red-black () 46 XYZ 各方向の隣りと斜めを同時に参照する27 点ステンシルなので,8 色必要 for(int i=0; i<nrow; i++){ SYMGS } カラーを回す外側ループ追加で 3 重化 for(int ic=0; i<8; ic++){ リカレンスがないのでディレクティブ 挿入で並列化 #pragma omp parallel for for(int i=st; i<=ed ++){ 8 z y x } } 演算順序を変更しカラー毎に処理することで 並列化が可能に

24 演算順序の変更に伴い, メモリ上での再配置も行い, メモリアクセスを連続化 47 1 行分の情報 27 double 27 int N^1 N N^1 N & N^2 N^1 N メッシュ点のグループ メモリ上配置 カラー 1 カラー カラー 3 カラー 4 ここで説明した 8 色カラーリングは問題 の構造的なメッシュ構造を前提とした プログラミングを禁止するレギュレーショ ンに抵触する. そのためもっと一般的なカラーリング 手法に変更. メッシュ点を複数まとめてグループを 構成する. そのグループ毎にカラーリングを実施. カラー内ではグループ間でリカレンス は発生しない. そこをスレッド並列に利用.

25 49 & ComputeSYMGS 144^147 OK ComputeSPMV 87^88 ComputeSYMGS 114^117 SPMV SPMV% MG% 20.0## 700.0## 18.0## 16.0## 14.0## 600.0## 500.0## Time%[s]% 12.0## 10.0## 8.0## 4.0 Time%[s]% 400.0## 300.0## ## 4.0## 2.0## 0.0## As#Is# Tune# 200.0## 100.0## 0.0## As#Is# Tune# ## Total % 700.0## 600.0## Time%[s]% 500.0## 400.0## 300.0## 200.0## 100.0## 0.0## As#Is# 19.2 Tune# ## ## ## 65536## 32768## 16384## GFLOPS'(log) 8192## 4096## 2048## 1024## GFLOPS# 512## 256## 128## 64## 32## 16## 8# 32# 128# 512# 2048# 8192# 32768# (log)

26 51 Titan Mira Site Computer Cores NSCC / Guangzhou RIKEN Advanced Inst for Comp Sci DOE/OS Oak Ridge Nat Lab DOE/OS Argonne Nat Lab Swiss CSCS Leibniz Rechenzentrum CEA/TGCC-GENCI Exploration and Production Eni S.p.A. DOE/OS L Berkeley Nat Lab Texas Advanced Computing Center Meteo France Meteo France U of Toulouse Cambridge U TiTech HPL Rmax (Pflops) HPL Rank HPCG (Pflops) HPCG/ HPL Tianhe-2 NUDT, Xeon 12C 2.2GHz + Intel Xeon 3,120, % Phi 57C + Custom K computer Fujitsu SPARC64 VIIIfx 8C + Custom Titan, Cray XK7 AMD 16C + Nvidia Kepler GPU 14C + Custom Mira BlueGene/Q, Power BQC 16C 1.60GHz + Custom Piz Daint, Cray XC30, Xeon 8C + Nvidia Kepler 14C + Custom 705, % 560, % 786, # 1.2% 115, % SuperMUC, Intel 8C + IB 147, % Curie tine nodes Bullx B510 Intel Xeon 8C 2.7 GHz + IB HPC2, Intel Xeon 10C 2.8 GHz + Nvidia Kepler 14C + IB Edison Cray XC30, Intel Xeon 12C 2.4GHz + Custom Stampede, Dell Intel (8c) + Intel Xeon Phi (61c) + IB Beaufix Bullx B710 Intel Xeon 12C 2.7 GHz + IB Prolix Bullx B710 Intel Xeon 2.7 GHz 12C + IB CALMIP Bullx DLC Intel Xeon 10C 2.8 GHz + IB Wilkes, Intel Xeon 6C 2.6 GHz + Nvidia Kepler 14C + IB TUSBAME-KFC Intel Xeon 6C 2.1 GHz + IB 79, % 62, % 132, # 2.7% 78, * % 24,192 23, (.467*).464 (.415*) % % 12, % % % HPL HPCG * scaled to reflect the sam number of cores # unoptimized implementa 52 ISC 14 では正式のベンチマークではなく参考のミニリストとして発表 京はこのミニリストで 2 位を獲得 1 位の天河 2 はピーク性能で京の 5 倍程度 HPCG の性能差は 1.3 倍 HPCG/HPL 性能比 :4.1% は群を抜いている HPCG は SC 14 で正式のベンチマークとなる予定

27 53 はじめに LINPACK HPCGとは HPCGベンチマークプログラム HPCGのチューニングと性能