http://www.mmsonline.com/articles/parallel-processing-speeds-toolpath-calculations TA : 菅 新 菅沼智史 水曜 新行紗弓 馬淵隼 木曜 情報知能工学演習V (前半第4週) 政田洋平 システム情報学研究科計算科学専攻
演習 V( 前半 ) の内容 第 1 週 : 高性能計算 (High Performance Computing = HPC) 向けプログラミングの基礎 第 2 週 : シミュレーションの基礎 第 3 週 : 波の移流方程式のシミュレーション 第 4,5 週 : OpenMP による並列化 第 6 週 : シミュレーションコードの OpenMP による並列化 到達目標 :HPC シミュレーション 並列化の基礎を理解する 授業予定 水曜 4/23, 4/30, 5/7, 5/14, 5/21, 5/28 木曜 4/24, 5/1, 5/8, 5/15, 5/22, 5/29, 6/5 参考図書 C/C++ プログラマーのための OpenMP 並列プログラミング菅原清文著
第 4 回の内容 先週の課題 ( 波の移流方程式 ) の解説並列化手法について (OpenMP, MPI, etc...) OpenMPの基礎 - forループの並列 - 共有 プライベート変数について
波の移流方程式ってどんな式?( 先週の復習 ) 移流方程式は波がその形を維持したまま発展していく方程式 : u t + c u x =0 u は波の流速 c は波面の速度 u u t = 0 t > 0 u2 u2 C u1 C u1 X X 移流方程式に従う波は... 初期条件としてある波の形 ( 空間分布 ) を与える. 厳密には波の形を保ったまま 波面が速度 C で移流する. ただし 数値的には波面が滑らかになっていく ( 鋭さがなくなる 鈍 [ なま ] る ).
先週の課題の解説 : 課題 1 ( 後退差分 v.s., 中心差分 ) 3 Wave Velocity U Wave Velocity U 2 1 0 789: -1-2 back central H#9: PH#9:Q@ARS;TUV1W 0 20 40 60 80 100 X Position X "#$%&'()*+,-./0123%456789:%"; 離散化手法の違いによって計算結果は変わる. この問題の場合 後退差分の方が式を安定に解く. どの離散化法が良いかは解きたい式の形による. ( 今日の資料の最後に付録 : フォンノイマンの安定性解析 )
先週の課題の解説 : 課題 2 ( 後退差分 v.s., 中心差分 ) 1 back20 back100 back1000 Wave Velocity U Wave Velocity U 0.8 0.6 0.4 0.2 exact N=20 N=100 N=1000 0 0 20 40 60 80 100 X Position X "#$%"&'()*+,-./0123456# 分割数 Nが大きくなるにつれ 計算精度は上昇. R#]^/0_`Fa)*3.$ST2/0%[\&6 良い結果を得るには 高解像度計算 が不可欠. "#$%"&'()*+,-./07839"6 分割数 Nが大きくなるにつれ 計算時間は増大. R#$ST2/0FU78VW%X)YZ*.OP/0%[\&6 高解像度 計算を 短時間で実行 したい. :;<=>:?@ABCDEFGHIJKLFMN-OPQ6# 並列化が必要
OpenMP による並列化 1
HPC でよく用いられる並列化手法 OpenMP (Open Multi-Processing) スレッド並列 ( 基本的に1 台の計算機内 ) ソースコードに指示行を挿入 MPI (Message Passing Interface) - 計算科学演習 (B4 or M1) - プロセス並列 ( 計算機間 計算機内どちらもOK) ソースコードにデータ通信の関数を追加数学ライブラリ並列化されたBLAS, LAPACK, ScaLAPACKなど行列演算 フーリエ変換などをサポート.
OpenMP (Open Multi-Processing) とは? 共有メモリ型計算機のためのプログラミング規格 C/C++ またはFortranのソースコードにディレクティブ ( 指示行 ) を挿入 C では #pragma omp のプラグマ指示行 shared memory multiprocessor CPU CPU CPU (thread) (thread) (thread) Memory cache cache cache bus
OpenMP の並列実行モデル Fork-join モデル #pragma omp parallel{ で囲まれた領域 = 並列リージョン 複数のスレッドで実行される並列実行を行うスレッドの集団 = チーム calculation1(); #pragma omp parallel { calculation2(); calculation3();
OpenMP のソースコード例 #pragma omp parallel を挿入するだけで その領域は並列化される. omp_get_thread_num() は 自分のスレッド番号 ( スレッド ID) を返す関数. #include <stdio.h> #include <omp.h> OpenMP int main(void){ #pragma omp parallel OpenMP { printf(" My thread number= %d \n", omp_get_thread_num());
コンパイル方法 C コンパイラにオプション fopenmp を付けるだけ gcc fopenmp ( ソースファイル名 ) 実行方法は普通の C プログラムと同じ./a.out #include <stdio.h> #include <omp.h> OpenMP int main(void){ #pragma omp parallel OpenMP { printf(" My thread number= %d \n", omp_get_thread_num()); オプション fopenmp をつけなければ #pragma 指示句は無視される ( 非 OpenMP プログラムとしてコンパイルされる )
スレッド数の制御 環境変数 OMP_NUM_THREADS で制御. スレッド数を変えて実行する場合は 環境変数を変えるだけ. 再コンパイルの必要はない. 8 スレッドの指定例 : tcsh の場合 : setenv OMP_NUM_THREADS 8 bash の場合 : export OMP_NUM_THREADS=8
練習問題 1: 次のソースコードをコンパイルして 1, 2, 4 スレッドで実行せよ. コンパイルオプションをつけない場合逐次実行されることを確認せよ. #include <stdio.h> #include <omp.h> int main(void){ #pragma omp parallel { printf(" My thread number= %d \n", omp_get_thread_num()); スレッド数は環境変数 OMP_NUM_THREADS で制御 コンパイルオプション -fopenmp
練習問題 2: 次のソースコードをコンパイルして 4 スレッドで実行せよ #include <stdio.h> #include <omp.h> int main(void){ #pragma omp parallel { printf(" First: My thread number= %d \n", omp_get_thread_num()); printf(" Second: My thread number= %d \n", omp_get_thread_num()); 2 つの printf 文は並列リージョンに入っているかいないかが違う. 並列リージョンに入っていない printf 文はマスタースレッド (thread_id = 0) が実行.
計算の並列化 ワークシェアリング (Work Sharing) 構文チーム内のスレッドで分担して実行する部分を指定並列リージョンで用いる for 構文ループを各スレッドで分担して実行 section 構文各セクションを各スレッドで分担して実行 master 構文マスタースレッドのみ実行
for 構文 (1) 並列リージョンの中で #pragma omp for を挿入する #include <stdio.h> #include <omp.h> int main(void){ int i; double x[100]; #pragma omp parallel { #pragma omp for for(i=0; i<100; i++){ x[i]=i; 直後の for ループを複数のスレッドで分割
for 構文 (2) ループ変数の分散について 例えば i = 0~~99 のループを 4 スレッドで分割する場合 均等分割 ( 最も効率よくなるように自動で分割 ) #pragma omp parallel { #pragma omp for for(i=0; i<100; i++){ x[i]=i;
練習問題 3: 1 右のソースコードをコンパイルして timeコマンドを使って 1,2,4スレッドでの計算時間を調べよ 2 別のターミナルを開いて topコマンドを実行した状態で 1,2,4スレッドの計算を行い CPU 使用率を調べよ #include <stdio.h> #include <omp.h> #include <math.h> int main(void){ int i, j; double x[1000000]; for(i=0; i<1000000;i++){ x[i]=0.0; for(j=0;j<250;j++){ #pragma omp parallel { #pragma omp for for(i=0; i<1000000; i++){ x[i]= x[i]+exp(i*0.00001)+exp(-i*0.0001); printf("%lf\n",x[999999]);
共有変数とプライベート変数 共有 (shared) 変数 : どのスレッドからも参照 更新が可能 プライベート (private) 変数 : それぞれのスレッドが異なる値を持つ変数 デフォルト値 : 指定をしない変数は基本的に共有変数 #pragma omp forの直後のループ変数はプライベート変数
共有 (shared) 変数 共有変数の指定方法は 並列化指示文の後に shared ( 変数名 ) を書く. 配列も共有される. 例 ) x,y,z を共有変数で宣言する場合 : #pragma omp parallel shared(x, y, z) 指定をしなければ 基本的に共有変数なので ほとんどの場合省略可.
プライベート (private) 変数 プライベート変数の指定方法は 並列化指示文の後 private ( 変数名 ) を書く. 配列もプライベート化される ( 各スレッドがそれぞれ配列を確保する ). for 文を並列化した場合 そのfor 文の変数は自動的に private 変数になる. 例 ) x,y,z をプライベート変数で宣言する場合 : #pragma omp parallel private(x, y, z) private 指定忘れが バグの原因になることが多い.
練習問題 4: #include <stdio.h> #include <omp.h> int main(void){ int thread_id; #pragma omp parallel private(thread_id) { thread_id= omp_get_thread_num(); sleep(1); // sleep (1 ) 次のコードの private(thread_id) を付けた場合と 付けない場合で値がどうなるか調べよ printf(" thread_id= %d, in %d\n", thread_id, omp_get_thread_num()); private 指定が無い場合は thread_id は共有変数になる. つまり全てのスレッドが上書き可能な変数になる.
ファーストプライベート (firstprivate) 変数 並列実行の開始時に逐次実行部分の変数の値を各スレッドの変数にコピーする. 並列実行中は private 変数として各スレッドがそれぞれの値を持つ. int n; n=1; #pragma omp parallel firstprivate(n) { printf("n=%d, Thread_ID=%d \n",n, omp_get_thread_num()); 並列リージョンに入るまでは共有変数 ( 共通の初期値を与える ). 並列リージョン突入後はプライベート変数として振る舞う.
練習問題 5: 次のコードの firstprivate(n) を付けた場合と private(n) を付けた場合で値がどうなるか調べよ #include <stdio.h> #include <omp.h> int main(void){ int n; n=1; #pragma omp parallel firstprivate(n) { printf("n=%d, Thread_ID=%d \n",n, omp_get_thread_num());
リダクション変数 並列計算時はそれぞれのスレッドで別々の値を持ち 並列リージョン終了時に各スレッドの値が足しあわさ れる ( 総和 ) 変数. 総和のほかに 積や最大値を求めることが可能. int i, sum; sum=0; #pragma omp parallel reduction(+:sum) { #pragma omp for for(i=0;i<10000;i++){ sum= sum + i; printf("sum=%d\n",sum);
練習問題 6: 次のベクトルの内積計算を行い その値を出力するコードを書き さらに for ループを並列化せよ ヒント : リダクション変数を使う
課題 4: 1 数値積分の並列化 : 次の f(x) について x=0.0 から 1.0 までの区間を積分計算し その値を出力するコードを作成せよ さらに for 文を OpenMP で並列化せよ 計算機上では 次式により近似値を計算する : 各長方形の幅は Δx xi はその中間地点 高さは f(xi) である 1,000,000 個の長方形の和として計算せよ exp の計算のためには #include <math.h> を付ければよい スレッド数を 1,2,4 で変化させた時 計算時間はどのように変化するか考察せよ. N = 1000000; for (i = 1; i <= N; i++){ ; ヒント : x = 0.5/(double)N + (double)i/(double)n ; sum = sum + exp(-x^2)/(double)n;
課題 4: x 2 + 10-5 sin(x)exp(x) 2 数値微分の並列化 : 次の f(x) について 区間 5.0 < x < 18.5 での微分値 df(x)/dx を計算し その絶対値の最大を出力するコードを作成せよ その際 for 文を OpenMP で並列化し 最大値は reduction 変数を用いて求めよ. f(x) = x 2 + 10-5 sin(x)exp(x) f(x) ただし 微分を求める各点の座標は x[i] = xleft + (xright - xleft)*(double)i/(double)n で与えられる. この時 分割数は N=10 3 xleft = 5.0, xright = 18.5 とせよ. 計算領域の左端と右端の微分値は計算に使用しない. 微分の離散化には中心差分を用いよ. xleft ヒント : 最大値の求め方の例 x xright
課題提出方法 ソースファイルを圧縮せずにメールに添付して提出メール本体には 学籍番号, 氏名を必ず明記締め切りは1 週間後の午後 1 時 ( 次回の授業の直前 ) 分からなかった点や授業の感想があれば書いてください提出先メールアドレスは ymasada_at_harbor.kobe-u.ac.jp [ _at_ @ ] Subject( 件名 ) は ENSHU5 :KADAI4: 学籍番号 ファイル名は ENSHU5 _KADAI4-1_ 学籍番号.c と ENSHU5 _KADAI4-2_ 学籍番号.c にはクラス (aまたはb), 学籍番号は半角小文字
( 補足資料 ) なぜ中心差分と後退差分で結果に違いが生じたのか?
計算アルゴリズムの数値的安定性 Von Neumann u n j = gn e ijθ :.. u n+1 i = u n i c t 2 x (un i+1 u n i 1) cδt/δx ν g = 1 1 2 ν(eiθ e iθ ) = 1 iνsinθ g 2 =1+ν 2 sin 2 θ 1 ν 1
計算アルゴリズムの数値的安定性 Von Neumann u n j = gn e ijθ : u n+1 i.. = u n i c t 2 x (un i u n i 1) cδt/δx ν g = 1 ν(1 e iθ ) g 2 = (1 ν + νcosθ) 2 + ν 2 sin 2 θ = 1 2ν(1 ν)(1 cosθ) 0 ν 1.