応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング

Presentation on theme: "応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング"— Presentation transcript:

1 応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング
応用数理工学特論 線形計算と ハイパフォーマンスコンピューティング 第2回 計算理工学専攻　張研究室 山本有作 日立製作所の山本有作です。 「～」について発表いたします。

2 前回の概要 １. マイクロプロセッサのメモリ階層 ２. 単体プロセッサ向けの高性能化手法 ３. アルゴリズムの工夫による演算量の削減
１.　マイクロプロセッサのメモリ階層 レジスタ，キャッシュ，主メモリ ２.　単体プロセッサ向けの高性能化手法 データ参照の局所性の向上 レジスタ有効利用のためのループ展開 キャッシュ有効利用のためのブロック化 行列乗算の例 ３.　アルゴリズムの工夫による演算量の削減 行列乗算のための Strassen アルゴリズム 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

3 今回の概要 １. 並列計算機の種類と特徴 ２. 共有メモリ型並列計算機 ３. 分散メモリ型並列計算機 プログラミングモデル 高性能化の技法
１.　並列計算機の種類と特徴 ２.　共有メモリ型並列計算機 プログラミングモデル 高性能化の技法 ３.　分散メモリ型並列計算機 本発表では，はじめに，研究の背景を述べてから，スパースソルバの概要，並列化手法，そして本研究で工夫した点の一つであるRISCプロセッサ向けの最適化についてご説明します。最後に，並列計算機SR2201上での性能評価とまとめを述べます。

4 １. 並列計算機の種類と特徴 共有メモリ型並列計算機 分散メモリ型並列計算機 SMPクラスタ 情報連携基盤センターの並列計算機
１.　並列計算機の種類と特徴 共有メモリ型並列計算機 分散メモリ型並列計算機 SMPクラスタ 情報連携基盤センターの並列計算機 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

5 スーパーコンピュータの性能動向 性能 年 2010 2000 1990 1980 1970 1960 1PFLOPS スカラー機
ASCI-5 性能 ベクトル機 ASCI-4 地球ｼﾐｭﾚｰﾀ ベクトル並列機 SR8000 並列機 VPP500 T3E-900 1TFLOPS CM-5 SR2201 nCUBE2 S3800 X-MP 1GFLOPS CRAY-1 S-810 スカラー／ベクトル機 → 並列機 CDC6600 (10 times faster / 20 years) 4 1MFLOPS IBM360/95 年 2010 2000 1990 1980 1970 1960

6 並列計算機の登場 並列計算機の普及の背景 並列計算機の特長 並列計算機の問題点 プロセッサの動作周波数向上の飽和
専用スーパーコンピュータの設計コストの増加 並列計算機の特長 プロセッサ数を増やすことでピーク性能を無制限に向上可能 分散メモリ型並列機では，プロセッサ数に比例した大きなメモリ空間が利用可能 汎用のプロセッサを使うことで設計コストを大幅に削減可能 並列計算機の問題点 多数のプロセッサを効率良く働かせるには，良い並列化アルゴリズムが必要

7 並列処理による高速化 複数のプロセッサで処理を分担することにより，プログラムの実行時間を短縮
プログラムの中で並列化対象部分の占める割合が大きいほど，高速化の効果が大きい １プロセッサ 並列計算機

8 さまざまな並列計算機 共有メモリ型並列計算機（SMP） SMPクラスタ 分散メモリ型並列計算機 地球シミュレータ Itaniumサーバ
Power Mac G5 分散メモリ型並列計算機 日立 SR11000 並列計算機の主流は，両者を融合させた形のSMPクラスタへ。 日立 SR2201 PCクラスタ Power Mac G5 クラスタ

9 チップレベルでの並列処理 最近では，ゲーム機の専用プロセッサが非常に高性能化 Xbox360用プロセッサ： 汎用PCの６倍の性能
PlayStation3用プロセッサ（Cell）：　汎用PCの20倍の性能 これらを数値計算に活用できれば，非常に低コストで 超高速の計算が可能 Cell ﾌﾟﾛｾｯｻ写真 Cell ﾌﾟﾛｾｯｻﾌﾞﾛｯｸ図 （9個のCPUを内蔵） PlayStation3

10 共有メモリ型並列計算機 構成 特徴 プログラミング言語 複数のプロセッサ（PU）がバスを 通してメモリを共有
　　通してメモリを共有 PUはそれぞれキャッシュを持つ。 特徴 メモリ空間が単一のためプログラミングが容易 PUの数が多すぎると，アクセス競合により性能が低下 　 →　４～８台程度の並列が多い。 プログラミング言語 OpenMP （FORTRAN/C/C++ ＋ 指示文）を使用 キャッシュ PU0 PU1 PU2 PU3 バス メモリ

11 分散メモリ型並列計算機 構成 特徴 プログラミング言語 各々がメモリを持つ複数のPUを ネットワークで接続 PUはそれぞれキャッシュを持つ。
　　ネットワークで接続 PUはそれぞれキャッシュを持つ。 特徴 数千～数万PU規模の並列が可能 PU間へのデータ分散を意識したプログラミングが必要 プログラミング言語 FORTRAN/C/C++ ＋ MPI を使用 キャッシュ PU0 PU1 PU2 PU3 メモリ ネットワーク

12 SMPクラスタ 構成 特徴 プログラミング 複数の共有メモリ型並列計算機 （SMP）をネットワークで接続
各ノードの性能を高くできるため，比較的少ないノード数で高性能を達成できる。 プログラミングは，ノード内部の計算では共有メモリ型並列機として，ノードをまたがる計算では分散メモリ型並列機として行う。 プログラミング MPI と OpenMP とを組み合わせて使用 PU0 PU1 PU2 PU3 メモリ PU0 PU1 PU2 PU3 メモリ PU0 PU1 PU2 PU3 メモリ ネットワーク

13 情報連携基盤センターの並列計算機 富士通 PrimePower HPC2500 （SMPクラスタ） 講義での利用 8GFLOPS／プロセッサ
64プロセッサ／ノード 全24ノード 全体で12TBの主記憶 講義での利用 「並列計算特論」の受講者にはアカウントを発行 受講が難しい人は山本まで連絡ください。 PU0 PU1 PU2 PU3 メモリ ネットワーク

14 並列アルゴリズムの例 (1) 数値積分によるπの計算 計算の並列化
π =∫0 1 4/(1+x2) dxの積分区間をn等分し，中点則により計算。 計算の並列化 n個の長方形を４個のプロセッサに割り当て，担当する長方形の面積の計算と部分和の計算を各プロセッサが行う。 各プロセッサからの寄与を合計する処理はプロセッサ０が行う。 4/(1+x2) x 1 プロセッサ0 プロセッサ1 プロセッサ2 プロセッサ3

15 並列アルゴリズムの例 (2) ２次元領域の温度変化の計算 計算の並列化 プロセッサ0 プロセッサ1
各格子点での温度変化を，隣接する４個の格子点との温度差から計算。 計算の並列化 格子を４個の領域に分割し，各領域に属する格子点での温度変化をその領域の担当プロセッサが計算。 プロセッサ2 プロセッサ3

16 分散メモリ型並列計算機での並列化 分散メモリ プロセッサ間通信 プロセッサ0 プロセッサ1 通信
各プロセッサは固有のメモリ空間を持ち，自分の担当する部分データを格納する。 共有メモリ方式に比べハードが作りやすく，超並列機に向く。 プロセッサ間通信 他プロセッサの持つデータを参照するには，通信が必要。 通信 プロセッサ2 プロセッサ3

17 分散メモリ型並列計算機での並列化（続き）
プロセッサ0 プロセッサ1 プログラム例 通信 PROGRAM HEAT REAL*8 A(4,4) ◆ 初期設定 DO ITER=1, 100 DO I=1, 4 DO J=1, 4 ◆ 必要なら隣接プロセッサ 　　　　　　　　よりAの値を受け取る ◆ A(I,J)の値を更新 END DO ◆ 結果の出力 STOP END プロセッサ2 プロセッサ3

18 並列化効率の向上 並列実行時間＝演算時間＋通信時間＋待ち時間 並列化効率＝ １プロセッサでの実行時間 pプロセッサでの実行時間 × p
プロセッサ０ プロセッサ１ プロセッサ２ プロセッサ３ 時間 ：　演算時間 ：　通信時間 ：　待ち時間

19 2.1 共有メモリ型並列計算機のプログラミング OpenMP とは OpenMP プログラムの構成要素 簡単な OpenMP プログラムの例
2.1　共有メモリ型並列計算機のプログラミング OpenMP とは OpenMP プログラムの構成要素 簡単な OpenMP プログラムの例 共有変数とプライベート変数 ループへのスレッド割り当ての指定 セクション型の並列化 参考文献 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

20 OpenMPとは 共有メモリ型並列計算機上での並列プログラミングのためのプログラミングモデル 米国のコンパイラメーカーを中心に仕様を決定
ベース言語（FORTRAN/C/C++）をディレクティブ（指示文）により並列プログラミングができるように拡張 米国のコンパイラメーカーを中心に仕様を決定 1997/10　FORTRAN Ver. 1.0 API 1997/10　C/C++ Ver. 1.0 API

21 OpenMPの実行モデル Fork-joinモデル 並列化を指定しない部分は逐次的に実行 指示文で指定された部分のみを複数のスレッドで実行
各スレッドは同じメモリ空間をアクセス 逐次実行部分 並列起動（fork） スレッド0 スレッド1 スレッド2 スレッド3 並列実行部分 並列終了（join） 逐次実行部分

22 OpenMP プログラムの構成要素 元のプログラム 指示文 ライブラリ関数 環境変数
通常の FORTRAN または C/C++ で書かれたプログラム 指示文 並列化すべき場所，並列化方法を指定 FORTRANでは !$OMP で始まる文 ライブラリ関数 並列実行部分でスレッド番号を取得するときなどに用いる。 環境変数 並列実行部分で使うスレッド数などを指定するのに用いる。 環境変数 OMP_NUM_THREADS でスレッド数を指定

23 簡単な OpenMP プログラムの例 (1) 環境変数 OMP_NUM_THREADS を 2 に設定 下記のプログラムをコンパイル・実行
実行結果 program hello integer omp_get_thread_num write(6,*) ‘program start.’ !$omp parallel write(6,*) ‘My thread number =’, omp_get_thread_num() !$omp end parallel write(6,*) ‘program end.’ Stop end スレッド番号を取得するライブラリ関数 指示文： 並列実行部分の開始 指示文：並列実行部分の終了 My thread number = 0 My thread number = 1

24 簡単な OpenMP プログラムの例 (2) 環境変数 OMP_NUM_HREADS を 2 に設定 下記のプログラムをコンパイル・実行
実行結果 1≦ i ≦50 がスレッド0，51≦ i ≦100 がスレッド1で計算される。 program axpy integer i double precision a, x(100), y(100), z(100) （a，x(100)，y(100)の値を設定） !$omp parallel do do i = 1, 100 z(i) = a*x(i) + y(i) end do （z(100)の値を表示） stop end 指示文： 直後の do ループの並列化指示 ベクトルの加算 z = ax + y

25 共有変数とプライベート変数 共有変数 プライベート変数
OpenMP のプログラミングモデルでは，基本的にすべての変数は共有変数（どのスレッドからも書き込み／読み出し可能） プログラム例(2)の a，x，y，z は共有変数の例 プライベート変数 ループインデックス変数 i については，スレッド 0 と 1 で共有すると，正しい制御ができない。そこで，各スレッドがそれぞれ別の変数を持つ必要がある。 このような変数をプライベート変数と呼ぶ。

26 変数の共有指定 デフォルト値 共有変数の指定（通常は不要） プライベート変数の指定 例(2)のプログラムで，指示を省略せずに書く場合
何も指示をしない変数については，基本的に共有変数となる。 しかし，並列化されたループのインデックス変数のように，明らかにプライベート変数でなければならない変数については，特に指示をしなくてもプライベート変数となる（ただし，多重ループの場合は並列化対象ループのインデックス変数のみ）。 共有変数の指定（通常は不要） 並列化指示文の後に，shared 節を追加する。 プライベート変数の指定 並列化指示文の後に，private 節を追加する。 例(2)のプログラムで，指示を省略せずに書く場合 !$omp parallel do shared(a, x, y, z) private(i)

27 変数の共有指定の例 ２重ループの並列化（行列ベクトル積） a, x, y は自動的に共有変数となる。 i は自動的にプライベート変数となる。
j はプライベート変数とすべきだが，自動的にはそうならない（並列化対象ループのインデックス変数ではない）ので指定が必要。 program gemv integer i, j double precision a(100,100), x(100), y(100) （a(100,100)，x(100)の値を設定） !$omp parallel do private(j) do i = 1, 100 y(i) = 0.0d0 do j = 1, 100 y(i) = y(i) + a(i,j) * x(j) end do （y(100)の値を表示） stop end j をプライベート変数に指定

28 リダクション変数 総和の計算（ベクトルの内積） 並列化の際，変数 c は共有変数とすべきか？ プライベート変数とすべきか？
そうすると，総和が正しく計算できない（排他制御の問題）。 プライベート変数とすべきか？ そうすると，並列実行部分終了後に，c の値が消えてしまう。 program dot integer i double precision x(100), y(100), c （x(100)，y(100)の値を設定） c = 0.0d0 do i = 1, 100 c = c + x(i)*y(i) end do （cの値を表示） stop end

29 リダクション変数（続き） リダクション変数とは 並列実行部分ではプライベート変数で，並列終了時に各スレッドの値が合計されるような変数
これにより，総和型のループも正しく並列化できる。 総和の他に，積や最大値を求めるリダクション変数指定も可能 program dot integer i double precision x(100), y(100), c （x(100)，y(100)の値を設定） c = 0.0d0 !$omp parallel do reduction(+:c) do i = 1, 100 c = c + x(i)*y(i) end do （cの値を表示） stop end c をリダクション変数に指定（総和型） ここでcの値が合計される。

30 ループへのスレッド割り当ての指定 指定なしの場合の割り当て方法 負荷の不均等が生じる場合
インデックス変数の値の範囲をスレッド数分に分割し，各部分をスレッドに割り当てる。 負荷の不均等が生じる場合 上三角行列とベクトルの積の計算などでは，この 　　やりかたで割り当てると負荷の不均等が発生 一定の長さLのブロックを周期的にスレッドに割り 　　当てるブロックサイクリック割り当てを採用すると， 　　負荷分散を改善できる場合がある。 × ＝ × ＝

31 ループへのスレッド割り当ての指定（続き）
ブロックサイクリック割り当ての例 上記の例は，コンパイル時に割り当てを決める静的割り当て 動的割り当てや，環境変数による割り当て方法指定も可能 program gemv integer i, j double precision a(100,100), x(100), y(100) （a(100,100)，x(100)の値を設定） !$omp parallel do private(j) schedule(static,10) do i = 1, 100 y(i) = 0.0d0 do j = i, 100 y(i) = a(i,j) * x(j) end do （y(100)の値を表示） stop end ブロック幅10のブロックサイクリック割り当て インデックス変数 j は陽にプライベート指定 上三角行列とベクトルの積

32 セクション型の並列化 今まで学んだのは，主に do ループの並列化
各スレッドに別々の仕事を並列に行わせたい場合は，セクション並列化を指定する。 program main （逐次処理の部分） !$omp parallel sections !$omp section （ここにスレッド0の処理内容を記述） （ここにスレッド1の処理内容を記述） !$OMP end parallel sections stop end セクション並列化の開始 セクション並列化の終了

33 参考文献 南里豪志，天野浩文： “OpenMP入門 (1), (2), (3)”，九州大学情報基盤センター広報，全国共同利用版, Vol.2, No.1 (2002), pp.1-40. 　　　/KOHO/OpenMP/intro.html 佐藤三久： “OpenMP”，JSPP’99チュートリアル, 1999. PGI Workstation User’s Guide

34 2.2 高性能化の技法 基本的な考え方 BLASの利用による高性能化
2.2　高性能化の技法 基本的な考え方 BLASの利用による高性能化 はじめに，研究の背景として，有限要素法とスパースソルバの必要性についてご説明します。

35 基本的な考え方 PU間の負荷分散均等化 PU間同期の削減 キャッシュメモリの有効利用 各PUの処理量が均等になるよう 処理を分割
　　処理を分割 PU間同期の削減 同期には通常，数百サイクルの時間が必要 並列粒度の増大による同期回数の削減が性能向上の鍵 キャッシュメモリの有効利用 演算順序の変更等により，キャッシュ中のデータの再利用性を向上 同時にPU間でのアクセス競合を削減し，メモリアクセスを高速化 キャッシュ PU0 PU1 PU2 PU3 バス メモリ 並列粒度とは，PU間での同期なしに行える処理の大きさ。

36 共有メモリ型並列計算機におけるメモリ階層の例
データ転送速度 非常に大 レジスタ 演算器 キャッシュ 8～128 ワード データ転送速度大 数100Kバイト ～ 数Mバイト データ転送速度小 主メモリ 数100Mバイト ～ 数Gバイト 演算器と主メモリの速度差は，年々増大。 キャッシュにデータがある間に演算をまとめて行う（データの再利用）ことにより，主メモリのデータ転送速度の遅さをカバーできる。

37 BLASの利用による高性能化 BLASとは BLASの種類 A A A C A B
Basic Linear Algebra Subprograms の略 個々のマシン向けにチューニングした基本行列演算のライブラリ BLASの種類 BLAS1:　ベクトルとベクトルの演算 内積　c := x t y AXPY演算　y: = ax + y など BLAS2:　行列とベクトルの演算 行列ベクトル積　y := Ax 行列のrank-1更新　A := A + xyt BLAS3:　行列と行列の演算 行列乗算　C := AB など A = × A A = ＋ × C A B = ×

38 BLASにおけるデータ再利用性と並列粒度
演算量：　O(N)，　データ量：　O(N) データ再利用性：　なし 並列粒度：　O(N/p)　（N： ベクトルの次元，p： プロセッサ台数） BLAS2 演算量：　O(N2)，　データ量：　O(N2) ベクトルデータのみ再利用性あり 並列粒度：　O(N2/p) BLAS3 演算量：　O(N3)，　データ量：　O(N2) O(N)回のデータ再利用が可能 並列粒度：　O(N3/p) 演算をできる限り BLAS2，3で行うことが高性能化のポイント A = × A A = ＋ × C B = A ×

39 現在利用可能な最適化BLAS プロセッサメーカーの提供するBLAS
Intel MKL，AMD ACML，IBM ESSL など メーカーによっては共有メモリ向け並列化版もあり ATLAS（Automatically Tuned Linear Algebra Subprograms） 対象プロセッサに合わせて自動的に性能を最適化するBLAS より入手可能 共有メモリ向け並列化版もあり GOTO BLAS テキサス大学オースチン校の後藤和茂氏によるBLAS より入手可能 この３種のBLASの中では，多くの場合最高速