第9回講義（2016年6月14日）-2

Transcription

1 スパコンプログラミング () (Ⅰ) 分解法 () 東京大学情報基盤センター准教授塙敏博 06 年 6 月 4 日 ( 火 )0:5-:0

2 06/6/7- スパコンプログラミング () (Ⅰ) 講義日程 ( 工学部共通科目 ). 4 月 9 日 ( 今日 ): ガイダンス. 4 月 6 日 並列数値処理の基本演算 ( 座学 ) 3. 5 月 0 日 : スパコン利用開始 ログイン作業 テストプログラム実行 4. 5 月 7 日 高性能プログラミング技法の基礎 ( 階層メモリ ループアンローリング ) 5. 5 月 4 日 高性能プログラミング技法の基礎 ( キャッシュブロック化 ) 6. 5 月 3 日 行列 - ベクトル積の並列化 レポートおよびコンテスト課題 ( 締切 : 06 年 8 月 8 日 ( 月 )4 時厳守 7. 6 月 7 日 (8:30-0:5) 大演習室 べき乗法の並列化 8. 6 月 7 日 (0:5-:0) 行列 - 行列積の並列化 () 9. 6 月 4 日 (8:30-0:5) 大演習室 行列 - 行列積の並列化 () 0. 6 月 4 日 (0:5-:0) 分解法 () コンテスト課題発表. 6 月 8 日 分解法 (). 7 月 5 日 分解法 (3) 3. 7 月 日 新しいスパコンの紹介 お試し 他

3 スパコンプログラミング () (Ⅰ) 3 分解法 ( 中級レベル以上 ) の演習日程 並列化が難しいので 3 週間確保してあります. 今週 講義 ( 知識 アルゴリズムの理解 ) 並列化の検討. 来週 分解法の逐次アルゴリズムの説明 分解法の並列化実習 () 3. 再来週 分解法の並列化実習 ()

4 スパコンプログラミング () (Ⅰ) 4 講義の流れ. 分解法 ガウス ジョルダン法 ガウス消去法 枢軸選択 分解法 外積形式 内積形式 クラウト法 ブロック形式ガウス法 縦ブロックガウ ス法 前進 後退代入. サンプルプログラムの実行 3. 並列化のヒント 4. 実習課題 5. レポート課題

5 スパコンプログラミング () (Ⅰ) 5 分解法の概略 いろいろな変種があります

6 6 密行列に対する連立一次方程式 以下の式 ここでは実数の密行列は実数のベクトルとすると 解ベクトル を求めること 解ベクトルを求める方法は 以下の二種類が知られている. 直接解法行列操作により厳密解を求める方法. 反復解法近似解を反復計算で解に収束させ求める方法

7 ガウス ジョルダン法 基本的な消去法により解を求める 第 ステップ 第 ステップ 最終ステップ 7 " 第一行をもとに係数を消去 " 第二行をもとに係数を消去 * * * * 割り算のみで解を得る

8 8 ガウス ジョルダン法 右辺 の代わりに単位行列 I を用意して同様の操作をすれば 最終ステップでは逆行列が求まる 各ステップでの計算量が同じなので 並列化時の負荷バランスが良い

9 ガウス消去法 対角線より上の要素をゼロにしない方法 第 ステップ 第 ステップ 最終ステップ 9 " 第一行をもとに係数を消去 " 第二行をもとに係数を消去この消去を前進消去 (forwrd eimitio) とよぶ * * " # #

10 ガウス消去法 前進消去後 最後の項から順に解を求めていく 0 この代入処理を 後退代入 (cwrd sustitutio) とよぶ * * " # # ) / ( / * *

11 ガウス消去法 ガウス消去法は ガウス ジョルダン法に比べ 消去演算をする範囲が少ない ( 基本行より下のみ ) 演算量が低下する : 3 ( / 3) 基本行より下のみ演算するため 並列化するとガウス ジョルダン法に比べて 負荷バランスの劣化を起こしやすい 並列処理に向かないと考えた専門家がいた 現在はデータ分散の改良や通信の隠蔽技法 ハードウエア能力向上から ガウス消去法のほうが高速である 3

12 ピボッティング ガウス ジョルダン法 ガウス消去法とも 基本行の係数がゼロだと ゼロによる除算が生じ 計算が続行できない これを回避するため 消去する列から最も係数の大きなものを選択して 基本行と入れ替える ( 枢軸選択 ピボッティング pivot seectio) " 第 行をもとに係数を消去 0

13 ピボッティング ピボッティングには以下の 種の方法がある. 完全ピボッティング更新対象全体から最大のものを選ぶ方法. 部分ピボッティング更新対象の列または行から最大のものを選ぶ方式 ピボッティングの手間 経験的な数値安定性から部分ピボッティングが用いられることが多い 3 "

14 分解法 ガウス消去法のような消去処理を行列演算として定式化 連立一次方程式の行列表記 : 4 " " " "

15 分解法 分解法では 以下の 3 つのステップで解を計算する 第 ステップ : 行列 の 分解 第 ステップ : 前進代入 第 3 ステップ : 後退代入 5 u u u u u " # # " ) ( ) ( c c c c c c " # : ベクトル c を求める c c c c u u u u u " # : 解ベクトル を求める

16 6 分解法 行列 の 分解 には データアクセスの違いから以下の3 種の方法が知られている. 外積形式ガウス法 (outer-product form) 普通の消去法から導出. 内積形式ガウス法 (ier-product form) 分解がなされたとして の対角要素を に固定して導出 3. クラウト法 (Crout method) 分解がなされたとして の対角要素を に固定して導出

17 7 分解法の種類 外積形式 (outer-product form) ガウス法 ガウス消去法と同等の操作で 分解する 第 列を消去したい場合 係数を用いて を消去

18 外積形式ガウス法 すなわち列の消去は これを行列表記にすると 行列 をとすると この消去は 8 i i i... ) / ( m " i i i... ) / (

19 外積形式ガウス法 一般的に したがって 分解は ここで はの要素の符号を反転させたものであり 容易に得られる 消去作業が終われば行列 が得られる 9 ) ( ) (

20 0 外積形式ガウス法 (C 言語 ) for (0; <; ) { } dtmp.0 / [][]; for (i; i<; i) { [i][] [i][]*dtmp; } for (j; j<; j) { dj [][j]; for (i; i<; i) { [i][j] [i][j] [i][]*dj; } 注意 : の対角要素は であることを仮定 ( 計算しない ) の対角要素を入れる 更新 参照

21 外積形式ガウス法 (Fortr 言語 ) do dtmp.0d0 / ( ) do i (i ) (i ) * dtmp eddo do j dj ( j) do i (i j) (i j) (i )*dj eddo eddo eddo 注意 : の対角要素は であることを仮定 ( 計算しない ) の対角要素を入れる 更新 参照

22 外積形式ガウス法のまとめ 外積形式ガウス法では分解列の右側の領域が更新される right-ooig アルゴリズムと呼ばれる 外積形式ガウス法は並列化に向く 処理の中心の更新領域が多い 負荷バランスよくデータ分散できる 更新処理が 分解行と分解列という少ないデータを所有するだけで 要素ごとに独立して行える

23 3 内積形式ガウス法 内積形式 (ier-product form) ガウス法 分解がなされたと仮定した上で 行列 の対角要素を として導出した方法 u # # # # # u が求まる が求まる u 3 u " 3 0 # u... 0 u u # " u u u

24 4 内積形式ガウス法 この導出作業を一般化すると 以下の二部分に分かれる (I) u の導出部 (II) (I) で得られた値を元に の導出部 まとめると (I) (II) i u u i ( i j i j ij u j ij u ) / j u ( i ( i 3... )... )

25 5 内積形式ガウス法 (C 言語 ) for (0; <; ) { } for (j0; j<; j) { } dj [j][]; for (ij; i<; i) { } [i][] [i][] [i][j]*dj; [][].0 / [][]; for (i; <; ) { [i][][i][]*[][]; } 参照 更新と参照更新

26 6 内積形式ガウス法 (Fortr 言語 ) do do j dj (j ) do ij (i ) (i ) (i j) * dj; eddo eddo ( ).0d0 / ( ) do i (i )(i ) * ( ) eddo eddo 参照 更新と参照更新

27 7 内積形式ガウス法のまとめ 内積形式ガウス法では 分解列の左側の領域が主に参照される eft-ooig アルゴリズムと呼ばれる 内積形式ガウス法の並列化 行列 を列方向分散 (*Cycic) 参照領域のデータがないので 通信多発 ( ベクトルリダクションが毎回必要 ) 行列 を行方向分散 (Cycic*) 上三角行列 の要素 ( データ数が少ない ) を所有すれば 独立して計算可能

28 8 クラウト法 クラウト法 (Cout Method) 分解がなされたと仮定した上で 行列 の対角要 素をとして導出した方法 (cf. 内積形式ガウス法 ) # u # u # 0 # # " 0 " # # u u u が求まる の第 列が求まる u

29 9 クラウト法 この計算を一般化すると の第 列を求める場合 i i j ij u j ( i... ) u の第 行を求める場合 j ( j i i u ij ) / ( j... )

30 30 クラウト法 (C 言語 ) [0][0].0/[0][0]; for (j; j<; j) { [0][j][0][j]*[0][0]; } for (0; <; ) { for (j0; j<; j) { dj[j][]; for (i; i<; i) { [i][][i][]-[i][j]*dj; } } [][].0/[][]; for (i0; i<; i) { di[][i]; for (j; j<; j) { [][j][][j]-di*[i][j]; } } for (j; j<; j) { [][j][][j]*[][]; } } 参照 参照更新更新

31 3 クラウト法 (Fortr 言語 ) ().0d0/() do j ( j) ( j) * ( ) eddo do do j dj(j ) do i (i )(i ) - (i j) * dj eddo; eddo ( ).0d0 / ( ) do i di( i) do j ( j)( j) di * (i j) eddo; eddo do j ( j)( j) * ( ) eddo eddo 参照 参照更新更新

32 3 クラウト法 クラウト法では 最内ループの交換ができる 長さ (-) のループ 長さ (-) のループの内 最も長いループを最内に移動可 ベクトル計算機で実行性能が良い 分解列および分解行の外側に つの参照領域 分散メモリ型並列計算機での実装が困難 どのようにデータ分割しても大量通信発生 共有メモリ型並列計算機では並列化可能 参照領域があれば分解列と分解行は独立に計算可能

33 ブロック形式ガウス法 行列 を小行列に分解し その小行列単位で 分解する方法 分解と行列 - 行列積で実現できる 具体的には ( 各小行列を各 が所有 ) とすると 右辺は

34 ブロック形式ガウス法 第 ステップ 第 ステップ 第 3 ステップ 分解 を転送 * を計算 を転送 * を計算 分解 を転送 3 を転送 を転送 3 を転送

35 35 ブロック形式ガウス法 対角要素が 分解して 行方向 列方向に部分的な 分解を転送する ブロック形式ガウス法の実現法は二通りある. 実際に小行列 の逆行列を求める方法例 ) -. 逆行列を求めず 分解を用いる方法例 ) の実装の場合 行列 - 行列積が主演算となる 高効率で実装可能

36 36 縦ブロックガウス法 縦ブロックガウス法は 列方向のみデータを分割する方法 (cf. ブロック形式ガウス法 ) 並列化した場合 内に列データを全て所有しているため ピボッティング処理が実装しやすい ブロック形式ガウス法は実装が難しい 外積形式ガウス法の並列化に比べ. 通信回数の削減. ループアンローリングによる性能向上が期待できる

37 37 縦ブロックガウス法 データアクセスパターン 参照 m- m- m- 更新 m- 並列更新 m- m-

38 38 縦ブロックガウス法 縦ブロックガウス法は ある幅ごとに 分解を行う この幅のことをブロック幅とよぶ ブロック幅を用いて設計されたアルゴリズムを一般的にブロック化アルゴリズムとよぶ ブロック化をすることで 演算カーネルが 重ループ ( レベル BS) から 3 重ループ ( レベル 3 BS) になる 実装による性能向上が得られやすい

39 39 縦ブロックガウス法 (C 言語 ) 実際のカーネル部分 for (jm; jm<m; jm) { for (jm; j<; j) { dj [jm][j]; for (ijm; i<; i) { [i][j][i][j] - [i][jm]*dj; } } } ループ jm j i についてループの展開 ( ループアンローリング ) 可能

40 40 縦ブロックガウス法 (C 言語 ) jm について 段のアンローリング for (jm; jm<m; m) { } for (jm; j<; j) { dj0 [jm ][j]; dj [jm][j]; for (ijm; i<; i) { [i][j][i][j] - [i][jm } } ]*dj0 - [i][jm]*dj;

41 4 縦ブロックガウス法 (C 言語 ) さらに j についても 段のアンローリング for (jm; jm<m; m) { for (jm; j<; j) { dj00 [jm ][j ]; dj0 [jm][j ]; dj0 [jm ][j]; dj [jm][j]; for (ijm; i<; i) { [i][j ][i][j ] -[i][jm ]*dj00 - [i][jm]*dj0; [i][j][i][j] -[i][jm ]*dj0 - [i][jm]*dj; } } } この処理は ループ内で 段 列分の消去を同時にしているとみなせる ( 多段多列同時消去法 )

42 4 縦ブロックガウス法 (Fortr 言語 ) 実際のカーネル部分 do jm m do jm dj (jm j) do ijm (i j) (i j) (i jm) * dj eddo eddo eddo ループ jm j i についてループの展開 ( ループアンローリング ) 可能

43 43 縦ブロックガウス法 (Fortr 言語 ) jm について 段のアンローリング do jm m- do jm dj0 (jm j) dj (jm j) & do ijm (i j) (i j) - (i jm ) * dj0 - (i jm) * dj eddo eddo eddo

44 44 縦ブロックガウス法 (Fortr 言語 ) さらに j についても 段のアンローリング do jm m- do jm dj00 (jm j ) dj0 (jm j ) dj0 (jm j) dj (jm j) do ijm (i j ) (i j ) - (i jm ) *dj00 & - (i jm) *dj0 (i j) (i j) - (i jm ) *dj0 & -(i jm) *dj eddo; eddo; eddo この処理は ループ内で 段 列分の消去を同時にしているとみなせる ( 多段多列同時消去法 )

45 45 縦ブロックガウス法 ブロック化するとできる通信隠蔽 縦ブロックガウス法において データを列方向ブロックサイクリック分散 (*Cycic(m)) するだけで実現可能 分解が必要なブロックを所有する. 優先して 分解を行い結果を放送. その他の行列更新を行う そのほかの. 分解データ受信待ち. 行列更新 通信と計算のオーバーラップ 通信時間隠蔽

46 代入計算 行列 を固定 右辺 を変えて計算する場合は前進代入 後退代入を並列化する必要がある 結論 : データ分散により 処理パターンは異なるが並列化可能 列方向分散方式 (*Boc) など ウエーブフロント処理で並列化 行方向分散方式 (Boc*) など 列単位で並列性 ( 放送処理が必要 )

47 スパコンプログラミング () (Ⅰ) 47 サンプルプログラムの実行 ( 分解法 )

48 スパコンプログラミング () (Ⅰ) 48 分解のサンプルプログラムの注意点 C 言語版 /Fortr 言語版のファイル名 -f.tr ジョブスクリプトファイル u.sh 中のキュー名を ecture から ecture8 ( 工学部共通科目 ) に変更し pjsu してください ecture : 実習時間外のキュー ecture8: 実習時間内のキュー

49 スパコンプログラミング () (Ⅰ) 49 分解法のサンプルプログラムの実行 以下のコマンドを実行する $ cp /home/z3005/-f.tr./ $ tr vf -f.tr $ cd 以下のどちらかを実行 $ cd C : C 言語を使う人 $ cd F : Fortr 言語を使う人 以下共通 $ me $ pjsu u.sh 実行が終了したら 以下を実行する $ ct u.sh.oxxxxxx

50 スパコンプログラミング () (Ⅰ) 50 分解法のサンプルプログラムの実行 (C 言語 ) 以下のような結果が見えれば成功 N 9 sove time [sec.] [MFOS] ss vue: e-07 Ccuted vue:.3057e-0 OK Test is pssed.

51 分解法のサンプルプログラムの実行 (Fortr 言語 ) 以下のような結果が見えれば成功 スパコンプログラミング () (Ⅰ) 5 NN 9 sove time[sec.] MFOS ss vue: Ccuted vue: OK Test is pssed.

52 スパコンプログラミング () (Ⅰ) 5 Fortr 言語のサンプルプログラムの注意 行列サイズ N( および プロセッサ数 NROCS) の宣言は 以下のファイルにあります u.ic 行列サイズ変数が NN となっています iteger NN prmeter (NN9)

53 スパコンプログラミング () (Ⅰ) 53 サンプルプログラムの説明 #defie N 9 数字を変更すると 行列サイズが変更できます #defie MTRIX 生成行列の種類の指定です にすると 枢軸選択なしでも解ける行列を設定します 以外 にすると 乱数で行列を設定します この行列を解くには 枢軸選択処理が必要です ( サンプルプログラムでは解けません ) 解の検査方法 解ベクトル が ベクトルとなるように の右辺 を計算して設定しています 残差ベクトルの ノルムが *N より大きくなるとエラーです

54 スパコンプログラミング () (Ⅰ) 54 サンプルプログラムの説明 MySove 関数の仕様 doue 型の密行列 と 右辺ベクトル を入力とします 分解を用いて を求解し 解ベクトル を出力します 分解のアルゴリズムは外積形式 (right-ooig) です その他 N9 の時の 分解後の行列 の値 およびベクトル c の値 (C 言語のもの ) が ファイル uc.dt にあります デバックに活用してください

55 スパコンプログラミング () (Ⅰ) 55 演習課題 MySove 関数を並列化してください 中級以上のレベルであり 簡単ではありません とりあえず N9 で並列化してください できたら N9 以上の大きな値にして実行してください N9 で動いても N384 で動かなくなることがあります これは おそらく 前進代入か 前進消去部分が間違っています 何が問題か分からなくなった時は. 分解後の の値を表示 OK なら. ベクトル c の値を表示 OK なら 3. ベクトル の値を表示 というように 段階を経て部分を特定し 地道にデバックしてください これは 並列プログラミングの鉄則です

56 並列化のヒント : データ分散方式 行列 およびベクトル c の計算担当領域は以下のようにすると簡単です ( それぞれ各 で重複して持ちます ) ( ただし以下は 4 の場合で 実習環境は 9 です ) 0 N/NROCS 3 対 通信関数 (MI_Sed MI_Recv) のみで実装できます 受信用バッファ (uf[n]) が必要です N N/ NROCS 0 3 スパコンプログラミング () (Ⅰ) 56 N/ NROCS c 0 3 N/ NROCS 0 3

57 スパコンプログラミング () (Ⅰ) 57 並列化のヒント : 分解部分 分解部分は 行列 に関して 最外の-ループがづつ変動し消去部分がづつ小さくなっていきます 現在のにおいて 対角要素から 行 ( 右図の青いベクトル 枢軸ベクトルと呼ぶ ) は 消去に必要な情報です 枢軸ベクトルなしでは 並列に消去できません 以上から 並列化する際 以下を考慮する必要があります. 対角要素を持っている 番号をどう計算するか. 対角要素を持っているは 担当範囲がつ小さくなる 3. 対角要素を持っているは 枢軸ベクトルを放送する ( その他のは受け取る ) N/NROCS 0 3 N

58 スパコンプログラミング () (Ⅰ) 58 並列化の道具 対角要素を持っている 番号は (*BOCK) 分割方式の場合で かつ - ループ ( 行目 ) の場合 以下のようになる. / i ここで i / umprocs; 枢軸ベクトルを放送する相手は 自分の 番号より大きく umprocs 番までの である

59 スパコンプログラミング () (Ⅰ) 59 並列化のヒント : 前進代入部分 前進代入部分は このデータ分散方法では 対角ブロック部分に相当するベクトルcの要素すべて決定し その後 対角ブロックに相当するベクトルcが各 で参照されます 対角ブロック部分の値が決定しないと 次の処理に進めませ ん N/ NROCS c 0 3 N/NROCS 0 3 N 0 3 N/ NROCS

60 スパコンプログラミング () (Ⅰ) 60 並列化のヒント : 前進代入部分 以上をまとめると :. 最外ループ は ブロック幅 i ごとに進みます. 対角ブロックを持っている は 対角ブロック用の計算 ( 注意 ) をして 対応する c の要素を確定します 対角ブロックを持っている の判定方法は 分解の場合と同じです 3. 対角ブロックをもつ は myid- から計算している c の部分を受け取り 計算後 myid に結果を送る 0 は受け取らない umprocs- は送らない 4. 対角ブロック担当 は 計算結果を送らない

61 スパコンプログラミング () (Ⅰ) 6 前進代入部分 : 処理の流れ ステップ N/ NROCS ステップ N/ NROCS c 0 3 c 0 3 確定 N/NROCS 送信 / 受信 N/NROCS N N N/ NROCS N/ NROCS

62 スパコンプログラミング () (Ⅰ) 6 前進代入部分 : 処理の流れ ステップ 3 N/ NROCS ステップ 4 N/ NROCS c 0 3 c 0 3 確定 N/NROCS N/NROCS 0 0 送信 受信 送信 3 3 N N N/ NROCS N/ NROCS

63 スパコンプログラミング () (Ⅰ) 63 後退代入部分 前進代入と同様な処理をします ただし後退代入は前進代入に比べ 以下の違いがあります. 後ろから処理が始まります. 対角ブロックでの 行列 の対角要素の割り算が必要です

64 スパコンプログラミング () (Ⅰ) 64 後退代入部分 ステップ N/ NROCS ステップ N/ NROCS 確定 N/NROCS N/NROCS 0 0 送信 / 受信 3 3 N N c 0 3 c 0 3 N/ NROCS N/ NROCS

65 スパコンプログラミング () (Ⅰ) 65 レポート課題. [0] MySove 関数を並列化せよ 各 で行列 について すべての範囲を確保してよい. [5] MySove 関数を並列化せよ 各 で行列 について 最低限の範囲を確保せよ 3. [30] MySove 関数を並列化せよ 枢軸選択処理を実装せよ 問題のレベルに関する記述 : 00: きわめて簡単な問題 0: ちょっと考えればわかる問題 0: 標準的な問題 30: 数時間程度必要とする問題 40: 数週間程度必要とする問題 複雑な実装を必要とする 50: 数か月程度必要とする問題 未解決問題を含む 40 以上は 論文を出版するに値する問題

66 スパコンプログラミング () (Ⅰ) 66 レポート課題 4. [30] MySove 関数を 同時多段多列消去法を用いて並列化せよ また 同時多段多列の個数 ( ブロック幅 ) をチューニングして 性能を評価せよ 5. [35] 4. に加え 各ループにアンローリングを施し 性能をチューニングせよ 6. [40] 5. に加え ノンブロッキング通信を用いて通信処理を高速化せよ. 分解 前進代入 後退代入処理において 通信と計算がオーバラップするようなアルゴリズムを採用せよ ここで前進代入 後退代入処理においては ウエーブフロント処理を考慮すること

67 スパコンプログラミング () (Ⅰ) 67 来週へつづく 分解法 ()