OpenMP/OpenACC によるマルチコア メニィコア並列プログラミング入門 Fortran 編第 Ⅱ 部 :OpenMP 中島研吾 東京大学情報基盤センター

Transcription

1 OpenMP/OpenACC によるマルチコア メニィコア並列プログラミング入門 Fortran 編第 Ⅱ 部 :OpenMP 中島研吾 東京大学情報基盤センター

2 2 OpenMP Login to Reedbush-U Parallel Version of the Code by OpenMP STREAM Data Dependency

3 3 Hybrid 並列プログラミング スレッド並列 + メッセージパッシング OpenMP+ MPI CUDA + MPI, OpenACC + MPI 個人的には自動並列化 +MPI のことを ハイブリッド とは呼んでほしくない 自動並列化に頼るのは危険である 東大センターでは現在自動並列化機能はコンパイラの要件にしていない ( 調達時に加点すらしない ) 利用者にももちろん推奨していない OpenMP が MPI より簡単ということはない データ依存性のない計算であれば, 機械的に OpenMP 指示文を入れれば良い NUMA になるとより複雑 :First Touch Data Placement

4 4 Flat MPI vs. Hybrid Flat-MPI:Each Core -> Independent core core core memory core core memory core core memory core core core core core Hybrid:Hierarchal Structure memory core core core core memory core core core core memory core core core core

5 5 Memory L1 L1 L1 L1 L1 L1 L1 L2 L1 L1 L1 L1 L1 L1 L1 L1 L1 C C C C C C C C C C C C C C C C HB M x N Number of OpenMP threads per a single MPI process Number of MPI process per a single node

6 6 並列プログラミングモデルによって各プロ セスの受け持つデータの量は変わる分散メッシュの数も各サイズも変わる example: 6 nodes, 96 cores Flat MPI pcube HB 4x pcube HB 16x pcube

7 OMP-1 7 共有メモリ型計算機 MEMORY C P U C P U C P U C P U C P U C P U C P U C P U SMP Symmetric Multi Processors 複数の CPU( コア ) で同じメモリ空間を共有するアーキテクチュア

8 OMP-1 8 OpenMP とは 共有メモリ型並列計算機用の Directive の統一規格 この考え方が出てきたのは MPI や HPF に比べると遅く 1996 年であるという 現在 Ver.4.0 背景 Cray と SGI の合併 ASCI 計画の開始 主な計算機ベンダーが集まって OpenMP ARB を結成し, 1997 年にはもう規格案ができていたそうである SC98 ではすでに OpenMP のチュートリアルがあったし, すでに SGI Origin2000 で OpenMP-MPI ハイブリッドのシミュレーションをやっている例もあった

9 OMP-1 9 OpenMP とは ( 続き ) OpenMP は Fortan 版と C/C++ 版の規格が全く別々に進められてきた Ver.2.5 で言語間の仕様を統一 Ver.4.0 では GPU,Intel-MIC 等 Co-Processor, Accelerator 環境での動作も考慮 OpenACC

10 OMP-1 10 基本的仕様 OpenMP の概要 プログラムを並列に実行するための動作をユーザーが明示 OpenMP 実行環境は, 依存関係, 衝突, デッドロック, 競合条件, 結果としてプログラムが誤った実行につながるような問題に関するチェックは要求されていない プログラムが正しく実行されるよう構成するのはユーザーの責任である 実行モデル fork-join 型並列モデル 当初はマスタスレッドと呼ばれる単一プログラムとして実行を開始し, PARALLEL, END PARALLEL ディレクティヴの対で並列構造を構成する 並列構造が現れるとマスタスレッドはスレッドのチームを生成し, そのチームのマスタとなる いわゆる 入れ子構造 も可能であるが, ここでは扱わない

11 OMP-1 11 Fork-Join 型並列モデル Master Master thread thread thread thread Master thread thread Master thread thread Master thread thread thread thread thread thread PARALLEL fork END PARALLEL join PARALLEL fork END PARALLEL join

12 OMP-1 12 スレッド数 環境変数 OMP_NUM_THREADS 値の変え方 bash(.bashrc) export OMP_NUM_THREADS=8 csh(.cshrc) setenv OMP_NUM_THREADS 8 たとえば,OMP_NUM_THREADS=4 とすると, 以下のように i=1~100 のループが 4 分割され, 同時に実行される do i= 1, 25 do i= 1,100 do i= 26, 50 do i= 51, 75 do i= 76, 100

13 OMP-1 13 OpenMP に関連する情報 OpenMP Architecture Review Board (ARB) 参考文献 Chandra, R. et al. Parallel Programming in OpenMP (Morgan Kaufmann) Quinn, M.J. Parallel Programming in C with MPI and OpenMP (McGrawHill) Mattson, T.G. et al. Patterns for Parallel Programming (Addison Wesley) 牛島 OpenMP による並列プログラミングと数値計算法 ( 丸善 ) Chapman, B. et al. Using OpenMP (MIT Press) 最新! 富士通他による翻訳 :(OpenMP 3.0) 必携!

14 OMP-1 14 OpenMP に関する国際会議 WOMPEI(International Workshop on OpenMP: Experiences and Implementations ) 日本 (1 年半に一回 ) WOMPAT( アメリカ ),EWOMP( 欧州 ) 2005 年からこれらが統合されて IWOMP となる, 毎年開催 International Workshop on OpenMP Eugene, Oregon, USA

15 OMP-1 15 OpenMP の特徴 ディレクティヴ ( 指示行 ) の形で利用 挿入直後のループが並列化される コンパイラがサポートしていなければ, コメントとみなされる

16 OMP-1 16 OpenMP/Directives Array Operations Simple Substitution!$omp parallel do do i= 1, NP W(i,1)= 0.d0 W(i,2)= 0.d0!$omp end parallel do Dot Products!$omp parallel do private(i)!$omp& reduction(+:rho) do i= is, ie RHO= RHO + W(i,R)*W(i,Z)!$omp end parallel do DAXPY!$omp parallel do do i= 1, NP Y(i)= ALPHA*X(i) + Y(i)!$omp end parallel do

17 OMP-1 17 OpenMP/Direceives Matrix/Vector Products!$omp parallel do private(i,j) do i= 1, N W(i,Q)= D(i)*W(i,P) do j= 1, INL(i) W(i,Q)= W(i,Q) + W(IAL(j,i),P) do j= 1, INU(i) W(i,Q)= W(i,Q) + W(IAU(j,i),P)!$omp end parallel do

18 OMP-1 18 OpenMP の特徴 ディレクティヴ ( 指示行 ) の形で利用 挿入直後のループが並列化される コンパイラがサポートしていなければ, コメントとみなされる 何も指定しなければ, 何もしない 自動並列化, 自動ベクトル化 とは異なる 下手なことをするとおかしな結果になる : ベクトル化と同じ データ分散等 (Ordering) は利用者の責任 共有メモリユニット内のプロセッサ数に応じて, Thread が立ち上がる Thread :MPI でいう プロセス に相当する 普通は Thread 数 = 共有メモリユニット内プロセッサ数, コア数 であるが最近のアーキテクチャでは Hyper Threading (HT) がサポートされているものが多い (1 コアで 2-4 スレッド )

19 OMP-1 19 メモリ競合 MEMORY C P U C P U C P U C P U C P U C P U C P U C P U 複雑な処理をしている場合, 複数のスレッドがメモリ上の同じアドレスにあるデータを同時に更新する可能性がある 複数の CPU が配列の同じ成分を更新しようとする メモリを複数のコアで共有しているためこのようなことが起こりうる 場合によっては答えが変わる

20 OMP-1 20 メモリ競合 ( 続き ) MEMORY C P U C P U C P U C P U C P U C P U C P U C P U 本演習で扱っている例は, そのようなことが生じないよう, 各スレッドが同時に同じ成分を更新するようなことはないようにする これはユーザーの責任でやること, である ただ多くのコア数 ( スレッド数 ) が増えるほど, メモリへの負担が増えて, 処理速度は低下する

21 OMP-1 21 OpenMP の特徴 ( 続き ) 基本は!omp parallel do ~!omp end parallel do 変数について, グローバルな変数と,Thread 内でローカルな private な変数に分けられる デフォルトは global 内積を求める場合は reduction を使う!$omp parallel do private(is,ie,i)!$omp& reduction(+:rho) do ip= 1, PEsmpTOT is= STACKmcG(ip-1) + 1 ie= STACKmcG(ip ) do i= is, ie RHO= RHO + W(i,R)*W(i,Z)!$omp end parallel do W(:,:),R,Z,PEsmpTOT などはグローバル変数

22 OMP-1 22 FORTRAN と C use omp_lib...!$omp parallel do shared(n,x,y) private(i) do i= 1, n x(i)= x(i) + y(i)!$ omp end parallel do #include <omp.h>... { #pragma omp parallel for default(none) shared(n,x,y) private(i) } for (i=0; i<n; i++) x[i] += y[i];

23 OMP-1 23 本講義における方針 OpenMP は多様な機能を持っているが, それらの全てを逐一教えることはしない 講演者も全てを把握, 理解しているわけではない 数値解析に必要な最低限の機能のみ学習する 具体的には, 講義で扱っている ICCG 法によるポアソン方程式ソルバーを動かすために必要な機能のみについて学習する それ以外の機能については, 自習, 質問のこと ( 全てに答えられるとは限らない )

24 OMP-1 24 最初にやること use omp_lib #include <omp.h> FORTRAN C 様々な環境変数, インタフェースの定義 (OpenMP3.0 以降でサポート )

25 OMP-1 25 OpenMP ディレクィヴ (FORTRAN) sentinel directive_name [clause[[,] clause] ] 大文字小文字は区別されない sentinel 接頭辞 FORTRAN では!$OMP, C$OMP, *$OMP, 但し自由ソース形式では!$OMP のみ 継続行には FORTRAN と同じルールが適用される 以下はいずれも!$OMP PARALLEL DO SHARED(A,B,C)!$OMP PARALLEL DO!$OMP+SHARED (A,B,C)!$OMP PARALLEL DO &!$OMP SHARED (A,B,C)

26 OMP-1 26 OpenMP ディレクィヴ (C) #pragma omp directive_name [clause[[,] clause] ] 継続行は \ 小文字を使用 ( 変数名以外 ) #pragma omp parallel for shared (a,b,c)

27 OMP-1 27 PARALLEL DO!$OMP PARALLEL DO[clause[[,] clause] ] (do_loop)!$omp END PARALLEL DO #pragma omp parallel for [clause[[,] clause] ] (for_loop) 多重スレッドによって実行される領域を定義し,DO ループの並列化を実施する 並び項目 (clause): よく利用するもの PRIVATE(list) SHARED(list) DEFAULT(PRIVATE SHARED NONE) REDUCTION({operation intrinsic}: list)

28 OMP-1 28 REDUCTION REDUCTION ({operator instinsic}: list) reduction ({operator instinsic}: list) MPI_REDUCE のようなものと思えばよい Operator +,*,-,.AND.,.OR.,.EQV.,.NEQV. Intrinsic MAX,MIN,IAND,IOR,IEQR

29 OMP-1 29 実例 A1: 簡単なループ!$OMP PARALLEL DO do i= 1, N B(i)= (A(i) + B(i)) * 0.50!$OMP END PARALLEL DO ループの繰り返し変数 ( ここでは i ) はデフォルトで private なので, 明示的に宣言は不要 END PARALLEL DO は省略可能 C 言語ではそもそも存在しない

30 OMP-1 30 実例 A2:REDUCTION!$OMP PARALLEL DO DEFAULT(PRIVATE) REDUCTION(+:A,B) do i= 1, N call WORK (Alocal, Blocal) A= A + Alocal B= B + Blocal!$OMP END PARALLEL DO END PARALLEL DO は省略可能

31 31 OpenMP 使用時に呼び出すことのできる 関数群 関数名 int omp_get_num_threads (void) int omp_get_thread_num (void) double omp_get_wtime (void) void omp_set_num_threads (int num_threads) call omp_set_num_threads (num_threads) 内容スレッド総数自スレッドのID MPI_Wtimeと同じスレッド数設定

32 OMP-1 32 OpenMP を適用するには?( 内積 ) VAL= 0.d0 do i= 1, N VAL= VAL + W(i,R) * W(i,Z)

33 OMP-1 33 OpenMP を適用するには?( 内積 ) VAL= 0.d0 do i= 1, N VAL= VAL + W(i,R) * W(i,Z) VAL= 0.d0!$OMP PARALLEL DO PRIVATE(i) REDUCTION(+:VAL) do i= 1, N VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO OpenMP ディレクティヴの挿入これでも並列計算は可能

34 OMP-1 34 OpenMP を適用するには?( 内積 ) VAL= 0.d0 do i= 1, N VAL= VAL + W(i,R) * W(i,Z) VAL= 0.d0!$OMP PARALLEL DO PRIVATE(i) REDUCTION(+:VAL) do i= 1, N VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO OpenMP ディレクティヴの挿入これでも並列計算は可能 VAL= 0.d0!$OMP PARALLEL DO PRIVATE(ip,i) REDUCTION(+:VAL) do ip= 1, PEsmpTOT do i= index(ip-1)+1, index(ip) VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO 多重ループの導入 PEsmpTOT: スレッド数あらかじめ INDEX(:) を用意しておくより確実に並列計算実施 ( 別に効率がよくなるわけでは無い )

35 OMP-1 35 OpenMP を適用するには?( 内積 ) VAL= 0.d0 do i= 1, N VAL= VAL + W(i,R) * W(i,Z) VAL= 0.d0!$OMP PARALLEL DO PRIVATE(i) REDUCTION(+:VAL) do i= 1, N VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO OpenMP ディレクティヴの挿入これでも並列計算は可能 VAL= 0.d0!$OMP PARALLEL DO PRIVATE(ip,i) REDUCTION(+:VAL) do ip= 1, PEsmpTOT do i= index(ip-1)+1, index(ip) VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO 多重ループの導入 PEsmpTOT: スレッド数あらかじめ INDEX(:) を用意しておくより確実に並列計算実施 PEsmpTOT 個のスレッドが立ち上がり, 並列に実行

36 OMP-1 36 OpenMP を適用するには?( 内積 ) VAL= 0.d0!$OMP PARALLEL DO PRIVATE(ip,i) REDUCTION(+:VAL) do ip= 1, PEsmpTOT do i= index(ip-1)+1, index(ip) VAL= VAL + W(i,R) * W(i,Z)!$OMP END PARALLEL DO 多重ループの導入 PEsmpTOT: スレッド数あらかじめ INDEX(:) を用意しておくより確実に並列計算実施 PEsmpTOT 個のスレッドが立ち上がり, 並列に実行 例えば,N=100,PEsmpTOT=4 とすると : 各要素が計算されるスレッドを指定できる INDEX(0)= 0 INDEX(1)= 25 INDEX(2)= 50 INDEX(3)= 75 INDEX(4)= 100

37 37 Matrix-Vector Multiply do i = 1, N VAL= D(i)*W(i,P) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*W(itemL(k),P) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*W(itemU(k),P) W(i,Q)= VAL

38 38 Matrix-Vector Multiply!$omp parallel do private(ip,i,val,k) do ip= 1, PEsmpTOT do i = INDEX(ip-1)+1, INDEX(ip) VAL= D(i)*W(i,P) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*W(itemL(k),P) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*W(itemU(k),P) W(i,Q)= VAL!$omp end parallel do

39 39 Matrix-Vector Multiply: Other Approach This is rather better for GPU and (very) many-core architectures: simpler structure of loops!$omp parallel do private(i,val,k) do i = 1, N VAL= D(i)*W(i,P) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*W(itemL(k),P) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*W(itemU(k),P) W(i,Q)= VAL!$omp end parallel do

40 40 omp parallel (do) omp parallel-omp end parallel はそのたびにスレッドを生成, 消滅させる :fork-join ループが連続するとこれがオーバーヘッドになることがある omp parallel + omp do/omp for!$omp parallel...!$omp do do i= 1, N...!$omp do do i= 1, N...!$omp end parallel 必須 #pragma omp parallel... #pragma omp for {... #pragma omp for {

41 41 OpenMP Login to Reedbush-U (separate file) [1] Parallel Version of the Code by OpenMP STREAM Data Dependency

42 42 OpenMP Login to Reedbush-U Parallel Version of the Code by OpenMP STREAM Data Dependency

43 43 ここでの目標 solve_pcg( 対角スケーリング, 点ヤコビ前処理付き CG 法 ) の並列化 solver_pcg.f/c (solve_pcg) OpenMP 指示行を入れるだけ (Index 使わない ) でやってみる

44 OMP-1 44 前処理付共役勾配法 Preconditioned Conjugate Gradient Method (PCG) Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end 実際にやるべき計算は : { } [ ] 1 z = M { r} 近似逆行列 の計算が必要 : [ ] 1 [ ] 1 M A, [ M ] [ A] 究極の前処理 : 本当の逆行列 [ ] 1 [ ] 1 M = A, [ M ] = [ A] 対角スケーリング : 簡単 = 弱い [ ] 1 [ ] 1 M = D, [ M ] = [ D]

45 OMP-1 45 対角スケーリング, 点ヤコビ前処理 前処理行列として, もとの行列の対角成分のみを取り出した行列を前処理行列 [M] とする 対角スケーリング, 点ヤコビ (point-jacobi) 前処理 [ M ] = D D solve [M]z (i-1) = r (i-1) という場合に逆行列を簡単に求めることができる 簡単な問題では収束する D 0 0 N D N

46 solve_pcg (1/3) do i= 1, N X(i) = 0.d0 W(i,2)= 0.0D0 W(i,3)= 0.0D0 W(i,DD)= 1.d0/D(i) ITR= N do L= 1, ITR!C!C !C {z}= [Minv]{r}!C !C=== do i= 1, N W(i,Z)= W(i,R)*W(i,DD)!C===!C!C !C RHO= {r}{z}!c !C=== RHO= 0.d0 do i= 1, N RHO= RHO + W(i,R)*W(i,Z)!C=== Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end

47 solve_pcg (2/3) 47!C!C !C {p} = {z} if ITER=1!C BETA= RHO / RHO1 otherwise!c !C=== if ( L.eq.1 ) then do i= 1, N W(i,P)= W(i,Z) else BETA= RHO / RHO1 do i= 1, N W(i,P)= W(i,Z) + BETA*W(i,P) endif!c===!c!c !C {q}= [A]{p}!C !C=== do i= 1, N VAL= D(i)*W(i,P) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*W(itemL(k),P) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*W(itemU(k),P) W(i,Q)= VAL!C=== Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end

48 solve_pcg (3/3) 48!C!C !C ALPHA= RHO / {p}{q}!c !C=== C1= 0.d0 do i= 1, N C1= C1 + W(i,P)*W(i,Q) ALPHA= RHO / C1!C===!C !C {x}= {x} + ALPHA*{p}!C {r}= {r} - ALPHA*{q}!C !C=== do i= 1, N X(i) = X(i) + ALPHA * W(i,P) W(i,R)= W(i,R) - ALPHA * W(i,Q) DNRM2= 0.d0 do i= 1, N DNRM2= DNRM2 + W(i,R)**2!C=== ERR = dsqrt(dnrm2/bnrm2) if (ERR.lt. EPS) then IER = 0 goto 900 else RHO1 = RHO endif r= b-[a]x DNRM2= r 2 BNRM2= b 2 Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end IER = 1 ERR= r / b

49 49 module solver_pcg contains solve_pcg (1/5) parallel computing by OpenMP subroutine solve_pcg & & ( N, NPL, NPU, indexl, iteml, indexu, itemu, D, B, X, & & AL, AU, EPS, ITR, IER, N2) use omp_lib implicit REAL*8 (A-H,O-Z) integer :: N, NL, NU, N2 real(kind=8), dimension(n) :: D, B, X real(kind=8), dimension(npl) :: AL real(kind=8), dimension(npu) :: AU integer, dimension(0:n) :: indexl, indexu integer, dimension(npl):: iteml integer, dimension(npu):: itemu real(kind=8), dimension(:,:), allocatable :: W integer, parameter :: R= 1 integer, parameter :: Z= 2 integer, parameter :: Q= 2 integer, parameter :: P= 3 integer, parameter :: DD= 4

50 solve_pcg (2/5) 50!$omp parallel do private(i) do i= 1, N X(i) = 0.d0 W(i,2)= 0.0D0 W(i,3)= 0.0D0 W(i,DD)= 1.d0/D(i)!$omp parallel do private(i,val,k) do i= 1, N VAL= D(i)*X(i) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*X(itemL(k)) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*X(itemU(k)) W(i,R)= B(i) - VAL BNRM2= 0.0D0!$omp parallel do private(i) reduction(+:bnrm2) do i= 1, N BNRM2 = BNRM2 + B(i) **2 Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end

51 solve_pcg (3/5) 51 ITR= N Stime= omp_get_wtime() do L= 1, ITR!$omp parallel do private(i) do i= 1, N W(i,Z)= W(i,R)*W(i,DD) RHO= 0.d0!$omp parallel do private(i) reduction(+:rho) do i= 1, N RHO= RHO + W(i,R)*W(i,Z) if ( L.eq.1 ) then!$omp parallel do private(i) do i= 1, N W(i,P)= W(i,Z) else BETA= RHO / RHO1!$omp parallel do private(i) do i= 1, N W(i,P)= W(i,Z) + BETA*W(i,P) endif Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end

52 solve_pcg (4/5) 52!$omp parallel do private(i,val,k) do i= 1, N VAL= D(i)*W(i,P) do k= indexl(i-1)+1, indexl(i) VAL= VAL + AL(k)*W(itemL(k),P) do k= indexu(i-1)+1, indexu(i) VAL= VAL + AU(k)*W(itemU(k),P) W(i,Q)= VAL C1= 0.d0!$omp parallel do private(i) reduction(+:c1) do i= 1, N C1= C1 + W(i,P)*W(i,Q) ALPHA= RHO / C1!$omp parallel do private(i) do i= 1, N X(i) = X(i) + ALPHA * W(i,P) W(i,R)= W(i,R) - ALPHA * W(i,Q) Compute r (0) = b-[a]x (0) for i= 1, 2, solve [M]z (i-1) = r (i-1) ρ i-1 = r (i-1) z (i-1) if i=1 p (1) = z (0) else β i-1 = ρ i-1 /ρ i-2 p (i) = z (i-1) + β i-1 p (i-1) endif q (i) = [A]p (i) α i = ρ i-1 /p (i) q (i) x (i) = x (i-1) + α i p (i) r (i) = r (i-1) - α i q (i) check convergence r end DNRM2= 0.d0!$omp parallel do private(ip,i) reduction(+:dnrm2) do i= 1, N DNRM2= DNRM2 + W(i,R)**2 ERR = dsqrt(dnrm2/bnrm2)...

53 solve_pcg (5/5) Stime = omp_get_wtime() do L= 1, ITR if (ERR.lt. EPS) then IER = 0 goto 900 else RHO1 = RHO endif Elapsed Time= Etime - Stime IER = continue Etime= omp_get_wtime() write (*,'(i5,2(1pe16.6))') L, ERR write (*,'(1pe16.6, a)') Etime-Stime, ' sec. (solver)' ITR= L deallocate (W) return end

54 OMP-1 54 >$ cdw ファイルコピー :Reedbush-U >$ cp /lustre/gt00/z30088/omp/omp-c.tar. >$ cp /lustre/gt00/z30088/omp/omp-f.tar. >$ tar xvf omp-c.tar >$ tar xvf omp-f.tar >$ cd multicore 以下のディレクトリが出来ていることを確認 omp stream >$ cd omp/src20 >$ make >$ cd../run >$ qsub go.sh

55 <$O-omp>/src20/Makefile parallel computing by OpenMP F90 = ifort F90OPTFLAGS= -O3 -qopenmp -ipo -xcore-avx2 -align array32byte F90FLAGS =$(F90OPTFLAGS).SUFFIXES:.SUFFIXES:.o.f.f90.c #.f90.o:; $(F90) -c $(F90FLAGS) $(F90OPTFLAG) $<.f.o:; $(F90) -c $(F90FLAGS) $(F90OPTFLAG) $< # OBJS = solver_pcg.o rcm.o struct.o pcg.o boundary_cell.o cell_metrics.o input.o main.o poi_gen.o pointer_init.o outucd.o TARGET =../run/sol20 all: $(TARGET) $(TARGET): $(OBJS) $(F90) $(F90FLAGS) -o $(TARGET) $(OBJS) $(F90FLAGS) clean: rm -f *.o $(TARGET) *.mod *~ PI* 55

56 56 ジョブ実行 実行方法 基本的にバッチジョブのみ インタラクティヴの実行は 基本的に できません 実行手順 ジョブスクリプトを書きます ジョブを投入します ジョブの状態を確認します 結果を確認します その他 実行時には 1 ノード (16 コア ) が占有されます 他のユーザーのジョブに使われることはありません

57 57 Reedbush-U ノードのブロック図 1 ノード :2CPU( ソケット ) 18 コア 基本的に 1 ソケット (1CPU) のみ使う 18 コア Socket #0 Socket #1

58 ジョブスクリプト (1/2)go1.sh /luster/gt00/t00xxx/multicore/omp/run/go1.sh スケジューラへの指令 + シェルスクリプト #!/bin/sh #PBS -q u-tutorial 実行キュー名 #PBS -N test ジョブ名称 ( 省略可 ) #PBS -l select=1:ncpus=18 ノード数, コア数 (1-18) #PBS -Wgroup_list=gt00 グループ名 ( 財布 ) #PBS -l walltime=00:05:00 実行時間 #PBS -e test.err エラー出力ファイル #PBS -o test.lst 標準出力ファイル C cd $PBS_O_WORKDIR 実行ディレクトリへ移動. /etc/profile.d/modules.sh 必須 58 export OMP_NUM_THREADS=18 スレッド数 (=ncpus, 1-18) export KMP_AFFINITY=granularity=fine,compact numactl./sol20 プログラム実行 numactl は必ず付けておく export KMP_AFFINITY=granularity=fine,compact 各スレッドがSocket#0の0 番から始まる各コアに順番に割り当てられる

59 ジョブスクリプト (1/2)go2.sh /luster/gt00/t00xxx/multicore/omp/run/go2.sh スケジューラへの指令 + シェルスクリプト #!/bin/sh #PBS -q u-tutorial 実行キュー名 #PBS -N test ジョブ名称 ( 省略可 ) #PBS -l select=1:ncpus=18 ノード数, コア数 (1-18) #PBS -Wgroup_list=gt00 グループ名 ( 財布 ) #PBS -l walltime=00:05:00 実行時間 #PBS -e test.err エラー出力ファイル #PBS -o test.lst 標準出力ファイル C cd $PBS_O_WORKDIR 実行ディレクトリへ移動. /etc/profile.d/modules.sh 必須 59 export OMP_NUM_THREADS=18 スレッド数 (=ncpus, 1-18) numactl./sol20 プログラム実行 numactl は必ず付けておく 各スレッドが Socket#0/#1 の各コアにランダムに割り当てられる

60 60 >$ cdw >$ cd multicore/omp/run >$ qsub go1.sh >$ cat test.lst ジョブ投入 INPUT.DAT NX NY NZ 1.00e e e-00 DX/DY/DZ 1.0e-08 EPSICCG

61 61 利用可能なキュー 以下の 2 種類のキューを利用可能 最大 8 ノードを使える u-lecture 8 ノード (288 コア ),10 分, アカウント有効期間中利用可能 全教育ユーザーで共有 u-tutorial 4 ノード (144 コア ),10 分, 講義 演習実施時間帯 lecture よりは多くのジョブを投入可能 ( 混み具合による ) 61 MPI Programming

62 バッチ処理とは スパコン環境では 通常は インタラクティブ実行 ( コマンドラインで実行すること ) はできません ジョブはバッチ処理で実行します バッチキュー ジョブの依頼 バッチ処理システムがジョブを取り出す 実行 ユーザ スパコン 62

63 バッチ処理を用いたジョブの実行方法 Reedbush システムにおいてバッチ処理は Altuir 社のバッチシステム PBS Professionul で管理されています ジョブの投入 : qsub< ジョブスクリプトファイル名 > #!/bin/bash #PBS -q u-lecture #PBS -Wgroup_list=gt00 #PBS -l select=8:mpiprocs=36 #PBS -l walltime=00:01:00 cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh mpirun./hello ジョブスクリプトファイルの例 キュー名 :u-lecture 利用グループ名 :gt00 63

64 バッチ処理システムの使い方 主要コマンド (Reedbush の場合 ) ジョブの投入 : qsub< ジョブスクリプトファイル名 > 自分が投入したジョブの状況確認 : rbstut 投入ジョブの削除 : qdel < ジョブ ID> バッチキューの状態を見る : rbstut--rsc バッチキューの詳細構成を見る : rbstut rsc-x 投げられているジョブ数を見る : rbstut-b 過去の投入履歴を見る : rbstut H 同時に投入できる数 / 実行できる数を見る :rbstut--limit 64

65 rbstat--rsc の実行画面例 $ rbstat --rsc QUEUE STATUS NODE u-debug [ENABLE,START] 54 u-short [ENABLE,START] 16 u-regular [ENABLE,START] ---- u-small [ENABLE,START] u-medium [ENABLE,START] u-large [ENABLE,START] u-x-large [ENABLE,START] 288 u-interactive [ENABLE,START] ---- u-interactive_1 [ENABLE,START] u-interactive_4 [ENABLE,START] 54 u-lecture [ENABLE,START] 54 u-lecture8 [DISABLE,START] 54 u-tutorial [ENABLE,START] 使えるキュー名 ( リソースグループ ) 現在利用可能か 利用可能ノード数

66 rbstat--rsc-x の実行画面例 $ rbstat --rsc -x QUEUE STATUS MIN_NODE MAX_NODE MAX_ELAPSE REMAIN_ELAPSE MEM(GB)/NODE PROJECT u-debug [ENABLE,START] :30:00 00:30:00 244GB pz0105,gcxx u-short [ENABLE,START] :00:00 02:00:00 244GB pz0105,gcxx u-regular [ENABLE,START] ---- u-small [ENABLE,START] :00:00 12:00:00 244GB gcxx,pz u-medium [ENABLE,START] :00:00 12:00:00 244GB gcxx ---- u-large [ENABLE,START] :00:00 12:00:00 244GB gcxx ---- u-x-large [ENABLE,START] :00:00 06:00:00 244GB gcxx u-interactive [ENABLE,START] ---- u-interactive_1 [ENABLE,START] :15:00 00:15:00 244GB pz0105,gcxx ---- u-interactive_4 [ENABLE,START] :05:00 00:05:00 244GB pz0105,gcxx u-lecture [ENABLE,START] :10:00 00:10:00 244GB gt00,gtyy u-lecture8 [DISABLE,START] :10:00 00:10:00 244GB gtyy u-tutorial [ENABLE,START] :10:00 00:10:00 244GB gt00 使えるキュー名 ( リソースグループ ) 現在利用可能か ノードの実行情報 課金情報 ( 財布 ) 実習では 1 つのみ 66

67 rbstat--rsc-b の実行画面例 $ rbstat --rsc b QUEUE STATUS TOTAL RUNNING QUEUED HOLD BEGUN WAIT EXIT TRANSIT NODE u-debug [ENABLE,START] u-short [ENABLE,START] u-regular [ENABLE,START] ---- u-small [ENABLE,START] u-medium [ENABLE,START] u-large [ENABLE,START] u-x-large [ENABLE,START] u-interactive [ENABLE,START] ---- u-interactive_1 [ENABLE,START] u-interactive_4 [ENABLE,START] u-lecture [ENABLE,START] u-lecture8 [DISABLE,START] u-tutorial [ENABLE,START] 使えるキュー名 ( リソースグループ ) 現在使えるか ジョブの総数 実行しているジョブの数 待たされているジョブの数 ノードの利用可能数 67

68 68 PCG 計算時間 : Etime-Stime NX=NY=NZ=100,go1.sh 実行時間 : 安定しない :5 回くらい測定して一番速い時間を採用する :numactl をはずすと安定するが遅い Parallel Performance (%) IFLAG=0, compact Thre ad # sec Speed-up Large 大 Thread # Granularity 粒度 Problem Size/Thread Small 小

69 go1.sh Only cores on a single socket used 69 #!/bin/sh #PBS -q u-tutorial #PBS -N test01 #PBS -l select=1:ncpus=18 1,2,4,8,12,16,18 #PBS -Wgroup_list=gt00 #PBS -l walltime=00:10:00 #PBS -e test1.err #PBS -o test1.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export KMP_AFFINITY=granularity=fine,compact export OMP_NUM_THREADS=18 1,2,4,8,12,16,18 numactl./sol20

70 go2.sh cores are randomly selected from 2 sockets 70 #!/bin/sh #PBS -q u-tutorial #PBS -N test01 #PBS -l select=1:ncpus=18 1,2,4,8,12,16,18 #PBS -Wgroup_list=gt00 #PBS -l walltime=00:10:00 #PBS -e test2.err #PBS -o test2.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export KMP_AFFINITY=granularity=fine,compact export OMP_NUM_THREADS=18 1,2,4,8,12,16,18 numactl./sol20

71 Parallel Performance (%) Results: Parallel Performance NX=NY=NZ=100, IFLAG=0 Measurement: 5 times, best case go2.sh is better if Thread # is more than 8 based on IFLAG=0, compact with 1 thread IFLAG=0, compact IFLAG=0 go1.sh go2.sh Thread #

72 Results: Fluctuation Rate NX=NY=NZ=100, IFLAG=0 Measurement: 5 times, (Worst-Best)/Best Core Allocation of go2.sh is random 72 Fluctuation (%) IFLAG=0, compact IFLAG=0 go1.sh go2.sh Thread #

73 Reedbush-U 1-node: 2-CPU s/sockets cores 18 cores Socket #0 Socket #1 Each Node of Reedbush-U 2 Sockets (CPU s) of Intel Broadwell-EP Each socket has 18 cores Each core of a socket can access to the memory on the other socket : NUMA (Non-Uniform Memory Access) Utilization of the local memory is more efficient So far, only a single socket has been used Let s utilize both sockets

74 74 Exercises Effect of problem size (NX, NY, NZ) Effect of Thread # (OMP_NUM_THREADS: 1-18)

75 75 OpenMP Login to Reedbush-U Parallel Version of the Code by OpenMP STREAM Data Dependency

76 OMP-3 76 何故 18 倍にならないか? 18 スレッドがメモリにアクセスすると,1 スレッドの場合と比較して, スレッド当り ( コア当り ) メモリ性能は低下 飽和 疎行列は memory-bound なためその傾向がより顕著 疎行列計算の高速化 : 研究途上の課題

77 77 疎行列 密行列 do i= 1, N Y(i)= D(i)*X(i) do k= index(i-1)+1, index(i) kk= item(k) Y(i)= Y(i) + AMAT(k)*X(kk) do j= 1, N Y(j)= 0.d0 do i= 1, N Y(j)= Y(j) + A(i,j)*X(i) X in RHS 密行列 : 連続アクセス, キャッシュ有効利用 疎行列 : 連続性は保証されず, キャッシュを有効に活用できず より memory-bound

78 OMP-3 78 GeoFEM Benchmark ICCG 法の性能 ( 固体力学向け ) SR11K/J2 SR16K/M1 T2K FX10 京 Core #/Node Peak Performance (GFLOPS) STREAM Triad (GB/s) B/F GeoFEM (GFLOPS) % to Peak LLC/core (MB) 疎行列ソルバー :Memory-Bound

79 STREAM benchmark メモリバンド幅を測定するベンチマーク Copy: c(i)= a(i) Scale: c(i)= s*b(i) Add: c(i)= a(i) + b(i) Triad: c(i)= a(i) + s*b(i) Double precision appears to have 16 digits of accuracy Assuming 8 bytes per DOUBLE PRECISION word Number of processors = 16 Array size = Offset = 0 The total memory requirement is MB ( 45.8MB/task) You are running each test 10 times -- The *best* time for each test is used *EXCLUDING* the first and last iterations Function Rate (MB/s) Avg time Min time Max time Copy: Scale: Add: Triad:

80 80 マイクロプロセッサの動向 CPU 性能, メモリバンド幅のギャップ

81 OMP-1 81 実行 :MPI バージョン >$ cdw >$ cd multicore/stream >$ mpiifort -O3 -xcore-avx2 -align array32byte stream.f o stream >$ qsub XXX.sh 18 cores 18 cores Socket #0 0 th -17 th cores Socket #1 18 th -35 th cores

82 82 #!/bin/sh s01.sh: Use 1 core (0 th ) #PBS -q u-tutorial #PBS -N stream #PBS -l select=1:mpiprocs=1 MPI Process #(1-36) #PBS -Wgroup_list=gt00 #PBS -l walltime=00:05:00 #PBS -e err #PBS -o t01.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export I_MPI_PIN_PROCESSOR_LIST=0 mpirun./impimap.sh./stream use 0th core

83 83 s16.sh: Use 16 cores (0-15 th ) #!/bin/sh #PBS -q u-tutorial #PBS -N stream #PBS -l select=1:mpiprocs=16 MPI Process #(1-36) #PBS -Wgroup_list=gt00 #PBS -l walltime=00:05:00 #PBS -e err #PBS -o t16.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export I_MPI_PIN_PROCESSOR_LIST=0-15 use 0-15th core mpirun./impimap.sh./stream

84 84 s32.sh: Use 32 cores (16 ea) #!/bin/sh #PBS -q u-tutorial #PBS -N stream #PBS -l select=1:mpiprocs=32 MPI Process #(1-36) #PBS -Wgroup_list=gt00 #PBS -l walltime=00:05:00 #PBS -e err #PBS -o t32.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export I_MPI_PIN_PROCESSOR_LIST=0-15,18-33 mpirun./impimap.sh./stream

85 85 #!/bin/sh s36.sh: Use 36 cores (ALL) #PBS -q u-tutorial #PBS -N stream #PBS -l select=1:mpiprocs=36 MPI Process #(1-36) #PBS -Wgroup_list=gt00 #PBS -l walltime=00:05:00 #PBS -e err #PBS -o t36.lst cd $PBS_O_WORKDIR. /etc/profile.d/modules.sh export I_MPI_PIN_PROCESSOR_LIST=0-35 mpirun./impimap.sh./stream

86 86 Results of Triad on a Single Node of Reedbush-U, Peak is GB/sec. Thread # GB/sec Speed-up

87 87 Exercises Running the code Try various number of threads (1-36) export I_MPI_PIN_PROCESSOR_LIST=0-7 export I_MPI_PIN_PROCESSOR_LIST=0-3,18-21 OpenMP-version and Single PE version are available Fortran,C

88 88 OpenMP Login to Reedbush-U Parallel Version of the Code by OpenMP STREAM Data Dependency

89 OMP-1 89 ICCG 法の並列化 内積 :OK DAXPY:OK 行列ベクトル積 :OK 前処理

90 OMP-1 90 前処理はどうするか? 対角スケーリングなら簡単 : でも遅い do i= 1, N W(i,Z)= W(i,R)*W(i,DD)!$omp parallel do private(i) do i = 1, N W(i,Z)= W(i,R)*W(i,DD)!$omp end parallel do!$omp parallel do private(ip,i) do ip= 1, PEsmpTOT do i = INDEX(ip-1)+1, INDEX(ip) W(i,Z)= W(i,R)*W(i,DD)!$omp end parallel do 64*64*64 METHOD= E E E-09 real 0m14.662s METHOD= E E E E E E-09 real 0m19.660s

91 OMP-1 91 前処理はどうするか? 不完全修正コレスキー分解 前進代入 do i= 1, N VAL= D(i) do k= indexl(i-1)+1, indexl(i) VAL= VAL - (AL(k)**2) * W(itemL(k),DD) W(i,DD)= 1.d0/VAL do i= 1, N WVAL= W(i,Z) do k= indexl(i-1)+1, indexl(i) WVAL= WVAL - AL(k) * W(itemL(k),Z) W(i,Z)= WVAL * W(i,DD)

92 OMP-1 92 データ依存性 : メモリの読み込みと書き 出しが同時に発生し, 並列化困難 不完全修正コレスキー分解 前進代入 do i= 1, N VAL= D(i) do k= indexl(i-1)+1, indexl(i) VAL= VAL - (AL(k)**2) * W(itemL(k),DD) W(i,DD)= 1.d0/VAL do i= 1, N WVAL= W(i,Z) do k= indexl(i-1)+1, indexl(i) WVAL= WVAL - AL(k) * W(itemL(k),Z) W(i,Z)= WVAL * W(i,DD)

93 OMP-1 93 前進代入 4 スレッドによる並列化を試みる do i= 1, N WVAL= W(i,Z) do k= indexl(i-1)+1, indexl(i) WVAL= WVAL - AL(k) * W(itemL(k),Z) W(i,Z)= WVAL * W(i,DD)

94 OMP-1 94 前進代入 4 スレッドによる並列化を試みる !$omp parallel do private (ip,i,k,val) do ip= 1, 4 do i= INDEX(ip-1)+1, INDEX(ip) WVAL= W(i,Z) do k= indexl(i-1)+1, indexl(i) WVAL= WVAL - AL(k) * W(itemL(k),Z) W(i,Z)= WVAL * W(i,DD)!$omp parallel INDEX(0)= 0 INDEX(1)= 4 INDEX(2)= 8 INDEX(3)=12 INDEX(4)=16 do i=1,4 do i=5,8 do i=9,12 do i=13,16 このような 4 スレッドが同時に実施される

95 OMP-1 95 データ依存性 : メモリへの書き出し, 読み込みが同時に発生 !$omp parallel do private (ip,i,k,val) do ip= 1, 4 do i= INDEX(ip-1)+1, INDEX(ip) WVAL= W(i,Z) do k= indexl(i-1)+1, indexl(i) WVAL= WVAL - AL(k) * W(itemL(k),Z) W(i,Z)= WVAL * W(i,DD)!$omp parallel INDEX(0)= 0 INDEX(1)= 4 INDEX(2)= 8 INDEX(3)=12 INDEX(4)=16 #0 thread #1 thread #2 thread #3 thread の部分にデータ依存性発生 (1 のときは下三角成分無し )

96 OMP-1 96 ICCG 法の並列化 内積 :OK DAXPY:OK 行列ベクトル積 :OK 前処理 : なんとかしなければならない 単純に OpenMP などの指示行 (directive) を挿入しただけでは 並列化 できない