カーネルベンチマークコード開発の目的エクサスケール規模のシミュレーションの核となる数値計算アルゴリズムの中で特に重要なものについて数値計算ライブラリ等を用いてそのコストを推定するためにカーネルベンチマークを作成し評価に使用する対象計算アルゴリズム固有値計算 ( 実数密行列標準固有値計

Size: px

Start display at page:

Download "カーネルベンチマークコード開発の目的エクサスケール規模のシミュレーションの核となる数値計算アルゴリズムの中で特に重要なものについて数値計算ライブラリ等を用いてそのコストを推定するためにカーネルベンチマークを作成し評価に使用する対象計算アルゴリズム固有値計算 ( 実数密行列標準固有値計"

さみらうなだ
5 years ago
Views:

1 カーネルベンチマークコードの開発について EigenExa について EigenExa ベンチマークコードについてベンチマーク結果に基づく性能推定について 2/24/2014 理化学研究所三上和徳 1

2 カーネルベンチマークコード開発の目的エクサスケール規模のシミュレーションの核となる数値計算アルゴリズムの中で特に重要なものについて数値計算ライブラリ等を用いてそのコストを推定するためにカーネルベンチマークを作成し評価に使用する対象計算アルゴリズム固有値計算 ( 実数密行列標準固有値計算 ) 3 次元 FFT 乱数生成器公開計画表カーネルベンチマークコード固有値計算 ( 実対称密標準 ) 公開中 3 次元 FFT 4 月公開見込み乱数生成 3 月公開見込み 2

固有値計算実数対称密行列の標準固有値計算 AICS で開発中の EigenExa をもとにベンチマークコードを作成 ( 公開版 version 2.1a) h;p://www.aics.riken.

3 固有値計算実数対称密行列の標準固有値計算 AICS で開発中の EigenExa をもとにベンチマークコードを作成 ( 公開版 version 2.1a) h;p:// 10 種類のテスト行列 (+ 利用者指定固有値 ) プロセスグループ形状の指定固有値のみ固有値 + 固有ベクトルなど ( ただし全モード計算 ) 3

4 3D-FFT FFTW 向けに紹介されている代表的な 3DFFT の測定コードをもとに各種 FFT 実装や分割方法による性能を推定する AICS で開発中の KMATH_FFT からベンチマークを作成中 (3 月末完成 ) 3D- FFT の主要なライブラリ FFTE, FFTW をカーネルとして開発中 3 次元空間の 1 軸 2 軸 3 軸分割の選択 4

5 乱数生成器メルセンヌツイスター乱数生成器を利用した乱数生成部分の性能評価ベンチマークコード MPI のグループでできるだけ重複のない乱数列を管理し内部状態管理の I/O 処理も含む AICS で開発中の KMATH_RANDOM のベンチマークコードを現在改良しカーネルベンチ提供用に調整中 5

EigenExa について理化学研究所ホームページ 2013 年 12 月 5 日 60 秒でわかるプレスリリースより京を使い世界最高速の固有値計算に成功行列の固有値計算では行列を簡単な形式 ( 形状 ) に変換しそれを中間形式として取り扱います理研の研究チームは帯行列 ( ゼロでない要素が対角線上に帯状に分布する行列 ) を中間形式に採用することによって

6 EigenExa について理化学研究所ホームページ 2013 年 12 月 5 日 60 秒でわかるプレスリリースより京を使い世界最高速の固有値計算に成功行列の固有値計算では行列を簡単な形式 ( 形状 ) に変換しそれを中間形式として取り扱います理研の研究チームは帯行列 ( ゼロでない要素が対角線上に帯状に分布する行列 ) を中間形式に採用することによって前処理の時間の削減を図った新しい計算アルゴリズムを考案しそれを基にした数学ソフト EigenExa( アイゲンエクサ ) を開発しました京の全プロセッサを用いて計算した結果世界最大規模の 100 万 100 万の行列の固有値計算が 1 時間以内で可能なこと確認しましたこれまでの地球シミュレーターの記録 (40 万 40 万の行列で 3 時間半 ) を大幅に上回りました京の高い計算能力と EigenExa の利用により数十万から 100 万程度の固有値を求める問題は 1 時間以内にできることが立証されました今後シミュレーションの規模を大幅に拡大することが可能になりますなお EigenExa はオープンソフトウエアとして公開され理研計算科学研究機構研究部門のホームページからダウンロードできます 6

7 固有値計算 - 密行列解法の位置づけゼロ要素を落とし, 大規模問題での使用メモリ量演算量を削減反復法が基本疎行列ベクトル積 (spmv) が性能を左右固有値計算密行列解法全固有値全固有対全体の数分の 1 のモード行列の全要素 (NxN) をゼロと考えずに扱う直接的解法メモリ使用量 O(N^2), 演算量 O(N^3) 疎行列解法超大規模問題少数固有モード特定区間モード最小最大モード疎行列解法の内部解法に密行列解法を使用高性能高品質な密行列向けソルバの必要性 7

8 EigenExa - 世界の競争相手 ELPA 1step vs 2step の議論 :1step が高速の場合が多い三重対角帯の変換部分は未だよい実装できず実装レベルでは困難か? B/Q でアセンブラチューニングの方向 &GPU 化へ DPLASMA PLASMA, MAGMA での 2step スキーム 1node はスケーラブルに動作 DAG タスクスケジューリング Eigen- Exa 新 1step スキーム採用京において 2^16 コアまでの動作確認通信コストの認識動作 & 通信モデル構築階層化アルゴリズム自動チューニングの取り込み GPU 版も準備へ ScaLAPACK ver 枠組みに大きな変化なし MPI が BLACS の標準にルーチンの強化 : ü 非対称行列ソルバー (PDHSEQR) ü 新ルーチン MRRR(PDSYEVR) 8

9 ScaLAPACK/DPLASMA テネシー大学 ICL h;p://icl.cs.utk.edu/dplasma/index.html ScaLAPACK : Scalable Linear Algebra PACKage a library of high- performance linear algebra rou^nes for parallel distributed memory machines. ScaLAPACK solves dense and banded linear systems, least squares, eigenvalue, and singular value problems. The key ideas includes a block cyclic data distribu^on for dense matrices and a block data distribu^on for banded matrices block- par^^oned algorithms to ensure high levels of data reuse well- designed low- level modular components DPLASMA : Distributed Parallel Linear Algebra So`ware for Mul^core Arch. DPLASMA is the leading implementa^on of a dense linear algebra package for distributed heterogeneous systems. It is designed to deliver sustained performance for distributed systems where each node featuring mul^ple sockets of mul^core processors, and if available, accelerators like GPUs or Intel Xeon Phi. DPLASMA achieves this objec^ve through the state of the art PaRSEC run^me, por^ng the Parallel Linear Algebra So`ware for Mul^core Architectures (PLASMA) algorithms to the distributed memory realm. 9

10 ELPA マックスプランク研究所 h;p://elpa- lib.di- berlin.mpg.de/wiki ELPA : Eigenvalue solvers for Petaflop Applica^ons Library ELPA is a Fortran- based high- performance computa^onal library for the (massively) parallel solu^on of symmetric or Hermi^an, standard or generalized eigenvalue problems. Once compiled, ELPA library rou^nes can be linked to from C, C++, Fortran etc. code alike. ELPA works as a "drop- in enhancement" for Scalapack- based infrastructures (arguably the de facto standard for high- performance parallel linear algebra). Thus, ELPA is not independent of this infrastructure, but rather builds on it. Necessary prerequisite libraries for ELPA (o`en already provided by HPC vendors) include: Basic linear algebra subrou^nes (BLAS) Lapack Basic linear algebra communica^on subrou^nes (BLACS) Scalapack 10

11 EigenExa - 国際競争力のある新規計算スキームの採用 tridiagonal eigenpairs ScaLAPACK DPLASMA 1step Scheme dense ELPA 2step Scheme (Byte/Flop が低い ) eigenpairs band 新 1step Scheme eigenpairs 高性能実装が困難全固有ベクトルを求める場合は 1step と大差ないという報告多数 Eigen- Exa 状況によっては逆変換 ( 三重対角帯 ) の高性能実装も視野にいれつつ, 良好なものを選択 11

12 固有値の計算パターン Ax = λx 行列のサイズが大きい場合方針 :A の固有値は相似変換 A - > P - 1 AP をしても不変簡単な P( 回転行列など ) で変換を多数回行い対角行列に収束させるより簡単な行列形式に中継的に変換して計算時間の短縮をはかる (1) 係数行列をより帯域の狭い中間行列に変換して (2) 中間行列の固有値を求め (3) 本来の係数行列の固有ベクトルを求め直す (1) の中間行列のパターン毎に手法がある EigenExa の標準版 ScaLAPACK : 三重対角行列 :Householder( 鏡像 ) EigenExa の高速版 : 帯行列 (narrow- band 法 ) Development of a High- Performance Eigensolver on a Peta- Scale Next- Genera^on Supercomputer System, Imamura etal, Progress in NUCLEAR SCIENCE and TECHNOLOGY, Vol. 2, pp (2011) 商用ソフトウエアでは Lanczos 法 : 三重対角がよく利用される行列のサイズが小さい場合 : Jacobi 法やベキ乗法などの簡便な方法で直接求めても可 12

EigenExa - Parallel performance: strong scalability [sec] 2000 N=50K 200

N=10K 20 Eigen- K(N=10K) ELPA2- development ELPA2- dev(n=10k)

13 EigenExa - Parallel performance: strong scalability [sec] 2000 N=50K 200 N=20K N=130K K computer@riken AICS OpenMP+MPI hybrid 8thread/1proc/1node N=10K 20 Eigen- K(N=10K) ELPA2- development ELPA2- dev(n=10k) ScaLAPACK(N=10K) Faster K 2K 4K 8K 16K 32K 64K Part of the results is obtained by using the K computer at the RIKEN Advanced Ins^tute for Computa^onal Science (Proposal number hp12017). 2013/11/20 13 [cores] 13

14 EigenExa ベンチマークコードについて Fortran で書かれた主プログラム実行プログラム名 eigenexa_benchmark 京 (FX10) X86 Intel X86 GNU BlueGeneQ 用の Makefile 有り make # libeigenexa.aのみ生成される make eigenexa_benchmark # ベンチマークコードを生成 mpirun コマンドなどで起動実行オプションをコマンドライン引数と入力ファイルで指定テスト用の係数行列を自動生成複数の行列 ( タイプサイズ ) を連続実行可能 EigenExa の求解ルーチン eigen_sx() 又は eigen_s () を呼び出して実行固有値計算に必要な計算資源を出力表示 ( テキスト ) 14

15 EigenExa ベンチマークコードについて実行プログラム名 eigenexa_benchmark 実行時のオプション (X86 Linux Intel 環境での例 ) $./eigenexa_benchmark -help eigenexa_benchmark [options] options: -h displays this help and exit -f input_file uses input_file default is./in -g mode sets the process grid as follows R, r MPI_COMM_WORLD row-major mode C, c MPI_COMM_WORLD column-major mode A, a MPI_COMM_SELF (embarrasingly parallel) 1, 2,... 9 splitted MPI_COMM_WORLD with the color=mod(rank,{number}) -x dimx dimy sets the cartecian shape (dimx, dimy) dimx <= dimy must be hold. 15

16 EigenExa ベンチマークコードについて行列や求解モードは input_file で指定する input_file のレコードフォーマットは以下!! Input file format!! N bx by mode matrix solver!! N : matrix dimension! bx : block width for the forward transformation! by : block width for the backward transformation! mode : solver mode { 0 : only eigenvalues }! { 1 : eigenvalues and corresponding eigenvectors}! { 2 : mode 1 + accuracy improvement for eigenvalues}! matrix : test matrix { 11 types, }! solver : { 0 : eigen_sx, new algorithm, faster on the K }! { 1 : eigen_s, conventional algorithm }!! if a line starts from '!', the line is treated as a comment!

17 EigenExa ベンチマークコードについてテスト用の係数行列 matrix パラメタ (0-9): 行列要素を自動生成以下の 10 種類の行列タイプから選択 matrix パラメタ (10): 行列要素をユーザが指定外部ファイル W.dat から読み込むベンチマークは input_file で複数の行列 ( タイプサイズ ) を指定して連続実行可能 Matrix type = 0 (Frank matrix) Matrix type = 1 (Toeplitz matrix) Matrix type = 2 (Random matrix) Matrix type = 3 (Frank matrix 2) Matrix type = 4 (W: 0, 1,..., n-1) Matrix type = 5 (W: sin(pai*5*i/(n-1)+eps^1/4)^3) Matrix type = 6 (W: MOD(i,5)+MOD(i,2)) Matrix type = 7 (W: same as Frank matrix) Matrix type = 8 (W: Uniform Distribution, [0,1)) Matrix type = 9 (W: Gauss Distribution, m=0,s=1) Matrix type = 10 (W: Read from the data file 'W.dat') 17

18 EigenExa ベンチマークコードについてテスト結果出力例 (X86 Linux Intel 環境での例 ) $ export OMP_NUM_THREADS=8 $ export I_MPI_FABRICS=shm:ofa $ mpirun -np 16 ${bin_path}/eigenexa-2.1/eigenexa_benchmark -x 4 4 INPUT FILE=IN ====================================================== Solver = eigen_sx / via penta-diagonal format Block width = 48 / 128 NUM.OF.PROCESS= 16 ( 4 4 ) NUM.OF.THREADS= 8 Matrix dimension = Matrix type = 0 (Frank matrix) Internally required memory = [Byte] mode 'X' :: mode 'A' + accuracy improvement Elapsed time = [sec] FLOP = Performance = [GFLOPS] 続く 18

19 EigenExa ベンチマークコードについてテスト結果出力例 (X86 Linux Intel 環境での例続き ) max w(i)-w(i).true / w.true = E *** Eigenvalue Relative Error *** : PASSED max w(i)-w(i).true = *** Eigenvalue Absolute Error *** : FAILED Do not mind it. Relative error is small enough A _{1}= epsilon= E-016 max Ax-wx _{1}/Ne A _{1}= E *** Residual Error Test *** : PASSED ZZ-I _{F}= E-012 *** Orthogonality Test *** : PASSED 19

20 EigenExa ベンチマークコードについてテスト結果出力例 (X86 Linux Intel 環境プロセス pin ステートも表示 ) $ export I_MPI_DEBUG=5 $ export I_MPI_FABRICS=shm:ofa $ mpirun -np 16 ${bin_path}/eigenexa-2.1/eigenexa_benchmark -x 4 4 [0] MPI startup(): Rank Pid Node name Pin cpu [0] MPI startup(): vsp25 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp25 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp27 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp27 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp29 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp29 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp10 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp10 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp11 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp11 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp12 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp12 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp16 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp16 {8,9,10,11,12,13,14,15} [0] MPI startup(): vsp17 {0,1,2,3,4,5,6,7} [0] MPI startup(): vsp17 {8,9,10,11,12,13,14,15}... OS( ジョブスケジューラ ) が捕捉したメモリ量 : Max Memory : 936 MB 20

21 EigenExa ベンチマークコードについて計算各フェイズでの通信待ち等さらに細かな情報を得たい場合は Makefile 中で以下のマクロを 1 に設定して make するとよい DEBUGFLAG = $(MACRO_D_PREFIX)TIMER_PRINT=0 実行時 stdout の例 calc (u,beta) mat-vec (Au) update (A-uv-vu) calc v E+000 v=v-(uv+vu)u UV post reduction E-002 COMM_STAT BCAST :: REDUCE :: REDIST :: E+000 GATHER :: E-002 TRD-BLK GFLOPS TRD-BLK-INFO

22 EigenExa ベンチマークコードについてテスト結果出力例 ( 京コンピュータでの例 ) $ export OMP_NUM_THREADS=8 $ mpiexec -n 16./eigenexa_benchmark -x 4 4 INPUT FILE=IN ====================================================== Solver = eigen_sx / via penta-diagonal format Block width = 48 / 128 NUM.OF.PROCESS= 16 ( 4 4 ) NUM.OF.THREADS= 8 Matrix dimension = Matrix type = 0 (Frank matrix) Internally required memory = [Byte] mode 'X' :: mode 'A' + accuracy improvement Elapsed time = [sec] FLOP = Performance = [GFLOPS] max w(i)-w(i).true / w.true = E *** Eigenvalue Relative Error *** : PASSED max w(i)-w(i).true = *** Eigenvalue Absolute Error *** : FAILED Do not mind it. Relative error is small enough. 以下カット 22

23 EigenExa ベンチマークコードについてベンダー統計ツールとの数値比較京コンピュータの場合 (fipp) 行列要素生成処理も含んだ表示 Elapsed(s) MFLOPS MFLOPS/PEAK(%) MIPS MIPS/PEAK(%) Application Process Process Process Process Process 5... similar 16 procs Mem throughput Mem throughput Elapsed(s) _chip(gb/s) /PEAK(%) SIMD(%) Application Process Process Process Process Process 5... similar 16 procs 23

24 ベンチマーク結果表示値に基づく性能推定について演算量 FLOP : EngenExa eigen_sx() の浮動小数点演算量 ~= C f x N 3 Matrix dimension = N メモリ量 Internally required memory :MPI プロセスあたりの合計配列サイズ ~= C m x N 2 ( 次ページ注参照 ) 計算時間 Elapsed ^me : 演算通信等の合計経過時間実効ピーク性能比 Performance [GFLOPS]: 固有値計算の性能 ( 全プロセス合計値 ) プロセス当たりの性能はデータレイアウトプロセスマッピングでほぼ決まるテスト結果を用いて推定 ( 同じプラットフォームで外挿 ) できることより大規模な行列を計算する場合の浮動小数点演算量メモリ量演算時間より多数のノード数を利用する場合の演算時間通信の待ち時間は実行条件への依存が大きく推定方法は一様ではない 24

利用上の留意点必要メモリ量の表示ベンチマークコードによるメモリ量の算出 (Internally required memory ) MPI プロセスあたりの最大メモリ量 ( ソースプログラムで静的に読める値 ) ジョブが実際に必要とするメモリ量はより多いテキストライブラリ MPI バッファなどのセグメントは入っていない OpenMP スレッド並列の場合はスレッド private

25 利用上の留意点必要メモリ量の表示ベンチマークコードによるメモリ量の算出 (Internally required memory ) MPI プロセスあたりの最大メモリ量 ( ソースプログラムで静的に読める値 ) ジョブが実際に必要とするメモリ量はより多いテキストライブラリ MPI バッファなどのセグメントは入っていない OpenMP スレッド並列の場合はスレッド private な変数領域が相当量必要ベンチマークコードによるメモリ量の表示 Internally required memory = [Byte] OSによる統計出力メモリ量 ( 京コンピュータfipp 表示の例 ) スレッド数 MAX MEMORY SIZE (USE) MiB (1.4GB) MiB (2.1GB) 計算結果計算の結果 a(n,n) は対角化されるが a(1,1) 成分は FLOP 値 ( 浮動小数点演算回数 ) a(2,1) 成分は経過時間の値で各々上書きされて返ってくるので注意 25

26 ダウンロードページ h;p:// 26

CPU Levels in the memory hierarchy Level 1 Level 2... Increasing distance from the CPU in access time Level n Size of the memory at each level 1: 2.2

CPU Levels in the memory hierarchy Level 1 Level 2... Increasing distance from the CPU in access time Level n Size of the memory at each level 1: 2.2 FFT 1 Fourier fast Fourier transform FFT FFT FFT 1 FFT FFT 2 Fourier 2.1 Fourier FFT Fourier discrete Fourier transform DFT DFT n 1 y k = j=0 x j ω jk n, 0 k n 1 (1) x j y k ω n = e 2πi/n i = 1 (1) n DFT