第三回 MPI 実践セミナー HPC システムズ株式会社新規事業企画室南本和秀 Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 1 rights reserved.

Size: px

Start display at page:

ともみこのえ
5 years ago
Views:

5 並列処理の簡単な例 Y=a[0]+a[1]+ +a[99] Σ 1CPU 4CPU プロセス0 プロセス1 プロセス2 プロセス3 a=a[0]+ +a[24] b=a[25]+ +a[49] c=a[50]+ +a[74] d=a[75]+ +a[99] Y=a+b+c+d 並列処理を行う部分 Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 5 rights

7 並列化方式の比較メリットデメリットマルチスレッド :OpenMP 共有メモリ方式メモリモデルが一般的コンパイラによる自動並列データの送受がメモリコピーの処理に置き換わる物理的な拡張性が低いプログラムが分散メモリ環境では動作しないマルチプロセス :MPI 分散メモリ方式大規模なシステム構成が可能実行ファイルが共有メモリ環境でも動作処理の分割法によってメモリモデルが複雑プログラミングが複雑プロセス同士の情報交換が多いと計算全体の性能が上がらない将来の理想象共有メモリ型分散メモリ型のハイブリッド開発は逐次で行いマルチプロセスで実行 2009/3/26 7 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

8 MPI とは (Message Passing Interface) 複数のプロセスがお互いにメッセージを送受信する仕組みを提供するための標準規格同一 ( 単一 ) のプログラムを複数のプロセスで同時に実行 (SPMD:Single Program Multi Data) 複数プロセス起動プロセス0 プロセス1 プロセス2 プロセス3 a=a[0]+ +a[24] b=a[25]+ +a[49] c=a[50]+ +a[74] d=a[75]+ +a[99] Y=a+b+c+d Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 8 プロセス間通信 rights

10 1 対 1 通信を行うための処理プロセス 0 プロセス 1 送信データ受信同一のプログラムを同時に実行データ RANK により処理を分岐受信送信 if(rank==0){ MPI_Send( sendbuff, 1, MPI_INT, 1, 1,MPI_COMM_WORLD); }else if(rank==1){ MPI_Recv( recvbuff, 1, MPI_INT, 0, 1,MPI_COMM_WORLD, status); } RANK1 へ sendbuffをsend RANK0 から recvbufferに Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 Recv 10 rights

11 送受信関数の概要 MPI_Send( sendbuff, 1, MPI_INT, 1, 1,MPI_COMM_WORLD); ( 送信バッファ, サイズ, 型, 宛先, タグ, コミュニケータ ) MPI_Recv( recvbuff, 1, MPI_INT, 0, 1,MPI_COMM_WORLD, status); ( 受信バッファ, サイズ, 型, 宛先, タグ, コミュニケータコミュニケータ, ステータス ) 型 :( 付録参照 ) MPI_INT, MPI_REAL 等送受信したいバッファの型により変えるタグ : 同じプロセス同士の通信であっても異なる処理を行うための識別子 2009/3/26 11 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

12 必須のインクルード mpi.h (Fortranの場合はmpif.h) 変数 int rank int nproc 処理 1MPI_Init 2MPI_Comm_size 3MPI_Comm_Rank 4MPI_Finalize MPI を使用するための必須事項自ランク番号総ランク数 MPI 環境の初期化プロセスの総数取得自ランク取得 MPI 終了処理 ( 例 ) RANK0 から RANK1 へ情報を送信するプログラム 1. #include <stdio.h> 2. #include mpi.h MPIのためのインクルードファイル int main(int argc, char* argv[]){ 5. int rank, nproc, status; 6. MPI_Status status; 7. int buffer, Message=25; 8. MPI_Init(&argc, &argv); 1 9. MPI_Comm_rank(MPI_COMM_WORLD, &rank);2 10. MPI_Comm_size(MPI_COMM_WORLD, &nproc); if(rank==0){ 14. MPI_Send( Message, 1, MPI_INT, 1, 1,MPI_COMM_WORLD); 15. }else if(rank==1){ 16. MPI_Recv(buffer, 1, MPI_INT, 0, 1,MPI_COMM_WORLD, &status) ; 17. } MPI_Finalize(); return 0; 22. } Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 12 rights

14 並列化プログラム 1. #include <stdio.h> 2. #include "mpi.h" 3. int main(int argc, char* argv[]){ 4. int arr[100]; 5. int itr; 6. int nproc; /* Number of proccess */ 7. int rank; /* Number of my rank */ 8. int buf[25]; 9. /* Initialize of MPI */ 10. MPI_Init(&argc, &argv); 11. MPI_Comm_rank( 12. MPI_COMM_WORLD, &rank); 13. MPI_Comm_size( 14. MPI_COMM_WORLD, &nproc); 15. /* 配列 arrにデータを格納 */ 16. for(itr=0;itr<100;itr++){ 17. arr[itr]=itr; 18. } 19. if(rank==0){ 20. /* 配列データを各ランクへ配信 */ 21. int dist; 22. for(dist=1; dist<4; dist++){ 25. } 52. int sum=0; 53. for(itr=0; itr<25; itr++){ 26. // bufに入った 25 個の整数型を 54. sum=sum+buf[itr]; 27. // dist 番のランクに向かって 55. } 28. // 特定のタグ1を付けて 29. // 送信 MPI_Send( buf, 25, MPI_INT, dist, 1, 31. MPI_COMM_WORLD); 32. } 33. for(itr=0;itr<25;itr++){ 34. buf[itr]=arr[itr]; 35. } 36. }else{ 37. /* 0 番以外のrank 番号を持つプロセスの処理 */ 38. // buf 配列に 39. // 25 個の 40. // int 型の配列を 41. // プロセス0 番から 42. // タグ1で 43. // 受け取る 44. MPI_Recv( buf, 25, MPI_INT, 0, 1, 45. MPI_COMM_WORLD, 46. MPI_STATUS_IGNORE); 47. } printf("rank=%d,sum=%d n",rank,sum); 57. // 各 RANKの結果を集計 // 58. if(rank == 0){ 59. int sum_other=0; 60. int src; 61. for(src=1; src<4; src++){ 62. MPI_Recv( 63. &sum_other, 1, MPI_INT, src, 2, 64. MPI_COMM_WORLD, 65. MPI_STATUS_IGNORE); 66. sum=sum+sum_other; 67. } 68. printf("summation=%d n",sum); 69. }else{ 70. MPI_Send(&sum, 1, MPI_INT, 0, 2, 71. MPI_COMM_WORLD); 72. } 73. MPI_Finalize(); return 0; 76. } 48. //!!!!!!!!!!!!!!!!!!!!!!!!!! // 49. // 並列計算される部分 // 23. for(itr=0;itr<25;itr++){ 50. //!!!!!!!!!!!!!!!!!!!!!!!!!! // 24. buf[itr]=arr[dist*25+itr]; 51. Copyright 2009/3/26 /* 配列の中身を足し合わせ (C) 2009 HPC SYSTEMS, */ Inc. All rights 14

15 基本的な集団通信集団通信関数を使用することで複数のプロセスがメッセージ交換を一斉に行う MPI_Scatter 1 つの送信元 RANK(Root) から全プロセスの受信バッファへ決まったサイズのデータを RANK が小さい順に格納 MPI_Gather 全 RANK の送信バッファから 1 つの宛先 RANK(Root) の受信バッファへ決まったサイズのデータを RANK が小さい順に格納送信バッファ RANK 0 RANK 0 受信バッファ RANK 0 RANK 1 RANK 2 RANK 3 受信バッファ送信バッファ Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 15 rights

16 集団通信を使った足し算プログラム 1. #include <stdio.h> 2. #include "mpi.h" 3. int main(int argc, char* argv[]){ 4. int arr[100]; 5. int itr; 6. int ansr[4]; 7. int sum=0; 8. int nproc; /* 全プロセス数 */ 9. int rank; /* 自分のランク番号 */ 10. int buf[25]; 11. MPI_Init(&argc, &argv); 12. MPI_Comm_rank(MPI_COMM_WORLD, &rank); 13. MPI_Comm_size(MPI_COMM_WORLD, &nproc); 14. if(rank==0){ 15. /* 配列 arrにデータを格納 */ 16. for(itr=0;itr<100;itr++){ 17. arr[itr]=itr; 18. } 19. } 24. sum=sum+buf[itr]; 25. } 26. printf("rank=%d,sum=%d n",rank,sum); 27. MPI_Gather(&sum, 1, MPI_INT, ansr, 1, MPI_INT, 0, 28. MPI_COMM_WORLD); 29. if(rank==0){ 30. sum=0; 31. for(itr=0; itr<4; itr++){ 32. sum=sum+ansr[itr]; 33. } 34. printf("sum=%d",sum); 35. } 36. MPI_Finalize(); 37. return 0; 38. } 20. MPI_Scatter(arr,25,MPI_INT, buf, 25, MPI_INT, 0, 21. MPI_COMM_WORLD); 22. /* 配列の中身を足し合わせ */ 23. for(itr=0; itr<25; itr++){ 2009/3/26 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights 16

17 MPI_Send/Recv に慣れよう! 一対一通信を使った具体例配列の総和を深く理解しよう 2009/3/26 HPC システムズ新規事業企画室 ( 工学博士 ) 渡邊啓正 2009 HPC SYSTEMS, Inc. All rights 1

18 お話しすること配列の総和 1/2 1/2 rank=0 1 1/nproc 1/nproc 1/nproc rank=0 1 nproc /3/ HPC SYSTEMS, Inc. All rights 2

19 rank=0 rank=1 main() { } 2009/3/26 MPI 初期化 Aの確保 Aの初期化 offset, numの計算 subaの確保 Aの後半をrank=1へ送信 Aの前半をsubAにコピー sumup(suba,num[0],&subtotal); answer = subtotal subtotal を rank=1 から受信 answer += subtotal 結果テスト MPI 終了処理メモリ解放 main() { } 2009 HPC SYSTEMS, Inc. All rights MPI 初期化 offset, numの計算 subaの確保 Aの後半をrank=0から受信 sumup(suba,num[1],&subtotal); subtotal を rank=0 へ送信 MPI 終了処理メモリ解放 3

20 nproc 並列版に変えるには (1/2) rank=1 から / への処理が nproc-1 倍に増える rank=0 rank=1..nproc-1 rank=1 への送信 rank=0 からの受信 rank=1..nproc-1への送信 ( 繰り返し ) rank=1 からの受信 rank=0 への送信 2009/3/26 rank=1..nproc-1からの受信 ( 繰り返し ) 2009 HPC SYSTEMS, Inc. All rights マスタワーカ型の場合 4

21 nproc 並列版に変えるには (2/2) C 言語でいうと rank=0 rank=1..nproc-1 for ( i = 1; i < nproc; i++ ) MPI_Send(A+offset[i], num[i], i ); MPI_Recv( A+offset[rank], num[rank], 0 ); for ( i = 1; i < nproc; i++ ) MPI_Recv(B+offset[i], num[i], i ); MPI_Send( B+offset[rank], num[rank], 0 ); 2009/3/ HPC SYSTEMS, Inc. All rights 5

rank=0 rank=1..nproc-1 main() { } 2009/3/26 MPI 初期化 Aの確保 Aの初期化 offset, numの計算 subaの確保 A の一部を rank=1..nproc-1 へ送信 A の一部を suba にコピー sumup(suba,num[0],&subtotal); answer = subtotal subtotalをrank=1.

22 rank=0 rank=1..nproc-1 main() { } 2009/3/26 MPI 初期化 Aの確保 Aの初期化 offset, numの計算 subaの確保 A の一部を rank=1..nproc-1 へ送信 A の一部を suba にコピー sumup(suba,num[0],&subtotal); answer = subtotal subtotalをrank=1..nproc-1 から受信 answer += subtotal 結果テスト MPI 終了処理メモリ解放 main() { } 2009 HPC SYSTEMS, Inc. All rights MPI 初期化 offset, numの計算 subaの確保 Aの一部をrank=0から受信 sumup(suba,num[rank],&subtotal); subtotal を rank=0 へ送信 MPI 終了処理メモリ解放 6

23 MPI プログラミング基礎編並列プログラムづくりのイメージを掴もう! 2009/3/26 HPC システムズ新規事業企画室 ( 工学博士 ) 渡邊啓正 2009 HPC SYSTEMS, Inc. All rights 1

24 お話しすることこうやれば MPI を使った並列プログラムを作れます! チューニング ( 効果的な高速化 )+ 並列化の技法知識と段取りの基本セット真似して作れるを目指す 2009/3/ HPC SYSTEMS, Inc. All rights 2

25 題材行列積 A B = C(matmul-seq.c) 105 for (row=0; row<n; row++) { for (col=0; col<n; col++) { sum = 0.0; for (k=0; k<n; k++) sum += A[row][k] * B[k][col]; C[row][col] = sum; } } N N N 複雑な実アプリでも並列化の基本は同じ A B = C 2009/3/ HPC SYSTEMS, Inc. All rights 3

26 並列化の一般的な流れ 1. ホットスポットの特定 2. 計算間の並列性の把握 3. データ分割方法粒度の検討 4. コード編集 5. コンパイルリンク 6. 実行 7. 性能評価 2009/3/ HPC SYSTEMS, Inc. All rights 4

27 1. ホットスポットの特定どの処理に一番時間がかかっているのか? そこの高速化がアプリ全体の高速化に効果的! 2009/3/26 main( ) {.;.; ;..;.; ;.; ;..; ; ;..; } かかった時間 HPC SYSTEMS, Inc. All rights 5

特定方法 1( 簡単 ): gprof を使う $ gcc pg foo.c $./a.out $ gprof a.out 関数ごとの実行時間統計 % cumulative self self total time seconds seconds calls ms/call ms/call name 77.78 0.07 0.07 4971903 0.00 0.00 mycompare 22.

28 特定方法 1( 簡単 ): gprof を使う $ gcc pg foo.c $./a.out $ gprof a.out 関数ごとの実行時間統計 % cumulative self self total time seconds seconds calls ms/call ms/call name mycompare myswap addtotable initresult mysort validateorder validatestability 何 % 占めるか実行時間何回呼ばれたか関数名 2009/3/ HPC SYSTEMS, Inc. All rights 6

29 特定方法 2( 難しい ): ソースコードから考えるソースコードを読んで演算量を見積もるプログラム部位ごとの int や float の演算回数 2009/3/26 プログラム部位演算回数 ( オーダー ) 入力受付 0 計算時間計測開始 0 メモリ確保 3 初期値設定 3N 2 行列積計算時間計測開始 0 計算時間表示 5 メモリ解放 0 8N 3 +3N HPC SYSTEMS, Inc. All rights N が大きいとき莫大な計算時間を要する 7

30 2. 計算間の並列性の把握並列性 : 複数の CPU へ分担させて同時に実行しても正しい計算結果が得られること行列積の場合 : 各 C[row][col] の計算は全て同時に実行して OK = = 同時実行可能! 2009/3/ HPC SYSTEMS, Inc. All rights = 8

31 並列性を考えた計算の実例領域分割計算熱伝導シミュレーション N 体問題粒子シミュレーション有限要素法の連立一次方程式解法パラメータスウィープ計算モンテカルロ計算気象予測金融リスク計算積分 ( 面積 ) の近似解生物データベース検索相互作用のない or 少ない領域に分けて計算を並列に実行する形多数の試行を並列に実行する形 2009/3/ HPC SYSTEMS, Inc. All rights 9

32 並列性で詰まったら ( その 1) データ依存性計算の順序に依って結果が変わる場合打開策 for (i=0;i<n;i++) X[i+1] = f(x[i]); 処理の大きさの単位並列の粒度を変えてみる処理の呼び出し元自体を並列化できないか? アルゴリズムを根本的に変える X[i+1] X[i+2] データ依存性が少ないアルゴリズムを誰かが作っているかもしれない論文調査など X[i] X[i+1] 2009/3/ HPC SYSTEMS, Inc. All rights 10

33 並列性で詰まったら ( その 2) フロー依存性前のフローに依って実行するかどうかが決まる場合打開策 a = BEFORE(); BEFORE if (a == 0) AFTER1(); else AFTER2(); 投機的実行 AFTER1 AFTER2 aの値に依らず AFTER1(); AFTER2(); を先走って実行させる aが決まり次第結果のどちらかを有効な値とする 2009/3/ HPC SYSTEMS, Inc. All rights 11

34 3. データ分割方法粒度の比較検討 C[row][col] の計算達を複数の CPU でどのように分担すればいいか時間のかかるネットワーク通信処理をできるだけ控えようローカルメモリアクセス時間 << ネットワーク通信時間データと計算作業をどのように分散させたらプロセス間のデータ通信回数通信量が最少データ通信回数通信量が最少になるか? = = 2009/3/ HPC SYSTEMS, Inc. All rights 12

35 分割 1 1x1 分割 ( 細切れ ) 通信回数が甚大です (4N 2 )!! データ通信に対して短い計算 = 非効率もっと大きめに分割 ( 分担 ) しましょうこりゃあかんわ = 2009/3/ HPC SYSTEMS, Inc. All rights 13

36 分割 2 A 横 B 縦ブロック分割 A を横に分割 B は縦に分割 C を横に分割分割はプロセス数分 (=P この図では 4) A[0] C C C C A[1] A[2] B [0] B [1] B [2] B [3] = C C C C C C C C A[3] C C C C 2009/3/ HPC SYSTEMS, Inc. All rights 14

37 分割 3 A 横ブロック分割 A を横に分割 B はコピー C を横に分割分割はプロセス数分 A[0] A[0] B=C[0] A[1] A[2] B = A[1] B=C[1] A[2] B=C[2] A[3] A[3] B=C[3] 2009/3/ HPC SYSTEMS, Inc. All rights 15

38 分割方法の比較分割方法通信回数通信量メモリ使用量 11x1 分割 4N 2 2N 3 +N 2 2N+1 2A 横 B 縦ブロック分割 2P 2-2 N 2 (P+3)(P-1)/P N 2 (2+1/P)/P 3A 横ブロック分割 3P-3 N 2 (P+2)(P-1)/P N 2 (1+2/P) P<<N 状況に応じてなにを重視するかで最適な戦略を考えましょう 2009/3/ HPC SYSTEMS, Inc. All rights 16

39 4. MPI コード製作 1. 主計算ルーチンの切り出し 2. MPIの定型句 3. 分割データ用メモリの確保 / 開放 4. データの分配 / 収集 5. 分割データを処理する主計算ルーチン 6. 並列計算結果のテストコード 7. 実行時間を計測するコード 2009/3/ HPC SYSTEMS, Inc. All rights 17

40 /3/26 main() { } (1) MPI 初期化 (2) A, B, Cの確保初期化 (3) 計測 1 (4) num[], offset[] の計算 (5) locala, localcの確保 (6) Aの分配 (7) Bのコピー (8) matmul_sub_parallel(locala, B, N/nproc, N, localc); (9) Cの収集 (10) locala, localcの解放 (11) num[], offset[] の解放 (12) 計測 2 (13) 計算結果のテスト (14) 2-1の表示 (15) A, B, Cの解放 (16) MPI 終了処理 2009 HPC SYSTEMS, Inc. All rights matmul-mpi.c 18

41 主計算ルーチンの切り出し 225 void matmul_sequential(const double *A, 226 const double *B, 227 const int size, 228 double *C) { 229 int row, col, k; 230 double sum; 231 for ( row = 0; row < size; row++ ) { 232 for ( col = 0; col < size; col++ ) { 233 sum = 0.0; 234 for ( k = 0; k < size; k++ ) { 235 sum += *(A + row * size + k) 236 * *(B + k * size + col); 237 } 238 *(C + row * size + col) = sum; 239 } 240 } 241 } 入出力を明確に後でテストに使う A B size matmul_ sequential C 2009/3/ HPC SYSTEMS, Inc. All rights 19

42 MPI の定型句を入れましょう (1/2) static const int root = 0; // Rank of the master process MPI_Init(&argc, &argv); // Initialize MPI system int nproc; // The number of processes being used MPI_Comm_size(MPI_COMM_WORLD, &nproc); // Calculate nproc int rank; // My process number MPI_Comm_rank(MPI_COMM_WORLD, &rank); // Calculate rank nproc=4 P rank=0(root) P rank=1 P rank=2 P rank=3 2009/3/ HPC SYSTEMS, Inc. All rights 20

43 MPI の定型句を入れましょう (2/2) char hostname[mpi_max_processor_name]; int namelen; MPI_Get_processor_name(hostname, &namelen); // Hostname fprintf(stdout, "[%d/%d] %s n", rank, nproc, hostname); fflush(stdout); // Output immediately ( 並列コード ) nproc=4 P rank=0 on node00 (CPU0) P rank=1 on node00 (CPU1) P rank=2 on node01 (CPU0) P rank=3 on node01 (CPU1) MPI_Finalize(); // Finalize MPI system 確認 : ちゃんと複数 CPU コアを使って動かせているか? 2009/3/ HPC SYSTEMS, Inc. All rights 21

44 分割データ用メモリの確保 / 開放 int *num = (int *)calloc(sizeof(int), nproc); int *offset = (int *)calloc(sizeof(int), nproc); for ( int i = 0; i < nproc; i++ ) { num[i] = (int)(n * N / nproc); offset[i] = (int)(n * N / nproc) * i; } 管理テーブル ( 今回は割り切れる場合のみ対応 ) N/nproc 0 N*N/nproc num[0] offset[0] double *locala = (double *)calloc(sizeof(double), num[rank]); double *localc = (double *)calloc(sizeof(double), num[rank]); N A[0] N*N free(locala); free(localc); free(num); free(offset); 各プロセスは AとC を横に nproc 分割した部分配列を持つ 2009/3/ HPC SYSTEMS, Inc. All rights 22

45 A の分配ルート (rank が 0 のプロセス ) が MPI_Send(A+offset[i], num[i], MPI_DOUBLE, i,...) を nproc 回繰返し, 各プロセスが次の受信を行う MPI_Recv(localA, num[rank], MPI_DOUBLE, root,...) MPI_Send x4 MPI_Recv rank A locala locala locala locala MPI_Recv MPI_Recv MPI_Recv A[0] 2009/3/ HPC SYSTEMS, Inc. All rights 23

46 B のコピールートが MPI_Send(B, N*N, MPI_DOUBLE, i,...) を nproc-1 回繰返し, 各プロセスが次の受信を行う MPI_Recv(B, N*N, MPI_DOUBLE, root,...) MPI_Send x3 rank B B B B B MPI_Recv MPI_Recv MPI_Recv 2009/3/ HPC SYSTEMS, Inc. All rights 24

47 C の収集全プロセスが MPI_Send(localC, num[rank], MPI_DOUBLE, root,...) を送信しルートが受信を nproc 回繰り返す MPI_Recv(C+offset[i], num[i], MPI_DOUBLE, i,...) MPI_Send rank localc localc localc localc MPI_Send MPI_Send MPI_Send A[0] B=C[0] C MPI_Recv x4 2009/3/ HPC SYSTEMS, Inc. All rights 25

48 分割データを処理する主計算ルーチン void matmul_sub_parallel(const double *A, const double *B, const int rowsize, const int size, double *C) { int row, col, k; 指定された行まで double sub; やれば OK for ( row = 0; row < rowsize; row++ ) { for ( col = 0; col < size; col++ ) { sub = 0.0; for ( k = 0; k < size; k++ ) { sub += *(A + row * size + k) * *(B + k * size + col); } *(C + row * size + col) = sub; } } } 2009/3/ HPC SYSTEMS, Inc. All rights 26

49 テストコード ( 逐次の計算結果との比較 ) double *testc = (double *)calloc(sizeof(double), N * N); matmul_sequential(a, B, N, testc); double epsiron = 1.0e-5; BOOLEAN isok = TRUE; int i; for (i = 0;i < N * N;i++) { double result = *(C + i); double test = *(testc + i); if (result > test) { if ((result - test) > epsiron) { isok = FALSE; break; } } else { if ((test - result) > epsiron) { isok = FALSE; break; } } 第三回 2009/3/26 free(testc); if (isok == FALSE) { printf("error: Result is NOT correct!! n"); } else { 2009 HPC SYSTEMS, Inc. All rights } printf("success: Result is correct!! n"); 結果が一つでも違っていたらバグあり 27

50 実行時間を計測するコード gettimeofday() で現在時刻を取得引き算 MPI_Barrier() を使ってプロセス全体をよーいどんさせてフェアに測定しよう rank 0 MPI_Barrier(MPI_COMM_WORLD); gettimeofday(&tv1, NULL); // Do something; MPI_Barrier(MPI_COMM_WORLD); gettimeofday(&tv2, NULL); double elapsed = 0.0; elapsed += (double)(tv2.tv_sec - tv1.tv_sec); elapsed += (double)((tv2.tv_usec - tv1.tv_usec) * 1e-6); /3/ HPC SYSTEMS, Inc. All rights 28

51 さらに深く知りたい人は実践 MPI-2 メッセージパッシングインタフェースの上級者向け機能ウイリアムグロップほか著株式会社ピアソンエデュケーション 2002 ISBN: 一通り MPI を使って並列プログラムを書けるようになった人向け入出力の最適化動的プロセス等 2009/3/ HPC SYSTEMS, Inc. All rights 29

52 5. コンパイル :mpicc > mpicc matmul-mpi.c MPI ライブラリが自動的にリンクされる gcc のオプション (-l や -o) をそのまま使える 2009/3/ HPC SYSTEMS, Inc. All rights 30

53 6. 実行 :mpirun > mpirun np 計算に使う CPU コア総数./a.out 並列計算の依頼ソフト LAN ユーザクラスタ結果複数 CPU を使って並列に計算 2009/3/ HPC SYSTEMS, Inc. All rights 31

54 7. 性能評価速度向上率並列プロセス数 vs. 速度向上 Linear matmul-mpi.c 多項式近似 /3/26 nproc top - 10:24:31 up 88 days, 12:11, 3 users, load average: 0.76, 0.26, 0.32 Tasks: 108 total, 2 running, 106 sleeping, 0 stopped, 0 zombie Cpu0 : 35.8% us, 27.5% sy, 0.0% ni, 36.7% id, 0.0% wa, 0.0% hi, 0.0% si Cpu1 : 43.9% us, 18.8% sy, 0.0% ni, 35.4% id, 0.0% wa, 0.4% hi, 1.3% si Cpu2 : 84.8% us, 15.2% sy, 0.0% ni, 0.0% id, 0.0% wa, 0.0% hi, 0.0% si Cpu3 : 33.9% us, 30.8% sy, 0.0% ni, 35.3% id, 0.0% wa, 0.0% hi, 0.0% si Mem: k total, k used, k free, 81608k buffers Swap: k total, 26192k used, k free, k cached PID USER PR NI VIRT RES SHR S %CPU %MEM TIME+ COMMAND hi-wata m 30m 2876 R :03.29 matmul-mpi.exe hi-wata m S :02.45 matmul-mpi.exe hi-wata m S :02.46 matmul-mpi.exe hi-wata m S :02.44 matmul-mpi.exe 2009 HPC SYSTEMS, Inc. All rights 32

55 休憩タイム 2009/3/26 第三回MPI実践セミナー実践セミナー Copyright (C) 第三回 2009 HPC SYSTEMS, Inc. All rights 33

56 MPI 環境の構築 MPI 環境を家で再現するには 2009/3/ HPC SYSTEMS, Inc. All rights 34

57 クラスタの機能プロセスと通信の制御例 : メッセージパッシングライブラリ LAN 接続セキュリティ情報の共有ファイルの共有クラスタの機能プロセスと通信の制御 LAN OS NIS NFS 2009/3/ HPC SYSTEMS, Inc. All rights 35

58 MPI のソフトウェア実装 OpenMPI メリット : CentOS4 で最初から導入済みきれいにコンポーネント分けされており, チューニングや障害解析がやりやすい複数ネットワークデバイスを自動認識する一斉に SSH で遠隔にアプリケーションを起動し, プロセスと通信路を管理する 2009/3/ HPC SYSTEMS, Inc. All rights 36

59 Platform Manager を使うと簡単! クラスタハードウェアに OS とミドルウェアを自動的にインストールしてくれるツールです簡単! オールインワンインタフェース 200 以上のワールドワイドユーザ企業対応 HW: Dell, HP, IBM, SGI, Sun 対応 OS: RHEL 3,4,5, SLES 9,10, CentOS App App ミドルウェア OS 2009/3/ HPC SYSTEMS, Inc. All rights 37

60 クラスタ構築の流れ 1 Linux PCにLinux DVD.isoファイルをダウンロード 2 Linux PCにPMをインストール 3 PMを起動 4 PM: クラスタイメージを作成 5 PM: クラスタイメージをサーバへ配置 6 クラスタを運用開始! 4 1 MPI 6 5 イメージ配置 LAN 2 PM 2009/3/ HPC SYSTEMS, Inc. All rights 38 3

61 GUI で簡単クラスタ構築ネットワーク構成 [New] ノード数クラスタ名詳しい資料は弊社スタッフまで 2009/3/ HPC SYSTEMS, Inc. All rights 39

62 VNC: リモートデスクトップツール VNC クライアントの起動遠隔に Linuxの GUI を使える 2009/3/ HPC SYSTEMS, Inc. All rights 40

63 コンパイル :mpicc > mpicc matmul-mpi.c MPI ライブラリが自動的にリンクされる gcc のオプション (-l や -o) をそのまま使えるホームディレクトリでコンパイル ( 推奨 ) 全ノードで共有されているため 2009/3/ HPC SYSTEMS, Inc. All rights 41

64 実行 :mpirun > mpirun --hostfile マシンファイルのパス np 計算に使う CPU コア総数./a.out マシンファイル OpenMPI で利用対象とするノードの情報 node01 slots=2 ノード上のプロセッサ数 node02 slots=2 node03 slots=2 node04 slots=2 /usr/local/openmpi1.2.8-intel91/etc/openmpi-default-hostfile 2009/3/ HPC SYSTEMS, Inc. All rights 42

65 おまけ :rank 割当て変更方法 > mpirun --hostfile マシンファイルのパス np 計算に使う CPU コア総数 --byslot./a.out --byslot 幅優先 --bynode 深さ優先 node01 node02 node03 node node01 node02 node03 node /3/ HPC SYSTEMS, Inc. All rights 43

66 画像ビューア :GQView S Slide show 2009/3/ HPC SYSTEMS, Inc. All rights 44

67 休憩タイム 2009/3/ HPC SYSTEMS, Inc. All rights 45

68 2 次元熱伝導熱伝導は基本的な物理現象の一つ熱伝導方程式から数値的に解析可能 T t = a 2 T 2 T ( x y ) a = λ /( ρc)[ m λ 熱伝導率 [ W ρ 密度 [ m c 比熱 [ J 3 ] 2 / s] / kgk] / mk] T 温度 [ K] t 時間 [ t] x, y 座標 [ m] ある場所とその近傍の温度を基にして少し未来の温度を知る Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 1 rights

69 差分解法差分解法差分解法差分解法 ) (, ) (, ) ( 1, ) ( 1, ) ( 1, ) ( 1, 1) (, 2 ) (, ) ( 1, ) ( 1, 2 ) (, ) ( 1, ) ( 1, ) (, 1) (, ) 4 ( ) 2 2 ( 2 / n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i n j i T T T T T T T y T T T x T T T a t T T h t h y x = = = = = γ γ とすると格子が等間隔であることから ) ( y T x T a t T + = 差分法で用いる形に変換差分法で用いる形に変換差分法で用いる形に変換差分法で用いる形に変換 2009/3/26 2 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

70 逐次版の動作 ( n+ 1) ( n) ( n) ( n) ( n) ( n) Ti, j γ ( Ti 1, j + Ti 1, j + Ti, j 1 + Ti, j 1 4Ti, j ) + T = + + T_new[i][j] = gamma*(te + Tw + Tn + Ts - 4.0*T[i][j]) + T[i][j]; Tw = T[i-1][j]; Te = T[i+1][j]; Ts = T[i][j-1]; Tn = T[i][j+1]; 本プログラムは大分大学工学部生産システム工学科三村泰成氏のアプリケーションを転載しました Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 3 rights

71 計算条件格子サイズは 256x256 境界以外の領域全体は 300 境界条件左端に幅 1 格子高さ 128 格子の熱源熱源の温度は 800 熱源以外の境界は断熱条件 [ 物質の性状 ] 密度 7.8e3; /* kg / (m^3) */ 熱伝導率 420.0; /* J/(kg K) */ 比熱 50.0; /* W/(m K) */ 2009/3/26 4 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

75 並列版のソースコード並列版のソースコード (pfdm2d_mpi.c 抜粋 ) 右の領域から左の領域へ情報を送る処理 if( rank!= 0 ) { /* send to left side */ for(j=0 ; j<nh ; j++) buf[j] = T[0][j]; MPI_Send(buf, NH, MPI_DOUBLE, rank-1, 1, MPI_COMM_WORLD); } if( rank!= nproc - 1 ) { /* receive from right side */ MPI_Recv(T_Right,NH,MPI_DOUBLE,rank+1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE); } 2009/3/26 8 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

76 穴埋め問題左の領域から右の領域へ情報を送る処理 if( rank!= nproc - 1) { /* send to right side */ for(j=0 ; j<nh ; j++) buf[j] = T[my_nw-1][j]; MPI_Send(, NH, MPI_DOUBLE,, 1, MPI_COMM_WORLD); } if(rank!= 0 ) { /* receive form left side */ MPI_Recv(,NH,MPI_DOUBLE,, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE); } 2009/3/26 9 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

77 アプリケーションの使い方解説実行方法 > mpirun np N./a.out M 4 コア 4 並列設定 :N=4, M=4 8 コア 8 並列設定 :N=8, M=8 time コマンドを併用すると全体の処理にかかる時間を計測することが可能 > time mpirun np N./a.out M 2009/3/26 10 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

78 < 実習 : コンパイル > > mpicc o pfdm2d pfdm2d_mpi.c alloc.c -o オプションは実行ファイルの名前を指定する (-o オプションなしにした場合デフォルト設定として同じディレクトリに a.out という実行ファイルが生成される ) > mpirun np 2 pfdm2d 2009/3/26 11 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

79 < 実習 : 答え合わせ > if( rank!= nproc - 1) { /* send to right side */ for(j=0 ; j<nh ; j++) buf[j] = T[my_nw-1][j]; MPI_Send(buf,NH, MPI_DOUBLE, rank+1, 1, MPI_COMM_WORLD); } if( rank!= 0 ) { /* receive form left side */ MPI_Recv(T_left, NH,MPI_DOUBLE, rank-1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE); } 2009/3/26 12 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

81 付録 1. MPI で用いる変数の型 ( 一部 ) MPI_CHAR MPI_BYTE MPI_SHORT MPI_INT MPI_LONG MPI_FLOAT MPI_DOUBLE MPI_UNSIGNED_CHAR MPI_UNSIGNED_SHORT MPI_UNSIGNED MPI_UNSIGNED_LONG char unsigned charと同様 short int long float double unsigned char unsigned short unsigned int unsigned long 2009/3/26 14 Copyright (C) 2009 HPC SYSTEMS, Inc. All rights

82 応用実習 2:N 体問題重力系シミュレーションを並列化しよう 2009/3/ HPC SYSTEMS, Inc. All rights 1

83 お話しすること N 体問題を MPI を使って並列化します! 基本アルゴリズム並列性の理解並列解法の理解 MPI_Send/Recvを書く ( 実習 ) 実行時間測定 ( 実習 ) 2009/3/ HPC SYSTEMS, Inc. All rights 2

84 N 体問題 N 個の質点間の相互作用の力の方程式を解いて行う物理シミュレーション一定時間後の速度と位置は? 計算量が莫大 ( 質点数の 2 乗 ) m m F= G r /3/26 銀河形成の研究等に研究等に重要 2009 HPC SYSTEMS, Inc. All rights 3

85 解説順序並列化すると? 2 体問題逐次版並列化前ベースを掴もう N 体問題逐次版 N 体問題並列版 2009/3/26 差分化された式って? 2009 HPC SYSTEMS, Inc. All rights 4

86 基本 :2 体問題のアルゴリズム質点 1 の万有引力の x 成分 : fx 2 ={ 万有引力の大きさ } {x の方向余弦 } = { G m1 m2 x1 x } { 2 r r 12 運動方程式を立てて離散化 ( 差分化 ): 質点 1の速度 : vk+ 1 = vk + fx2 / m1 * t x = x + v * t + v 質点 1 の位置 : 12 2 } k+ 1 k k k / 2* [x 1,y 1 ] m 1 [x 2,y 2 ] m 2 ここから質点が増えた数だけ万有引力が加算されると考えればよい t 2009/3/ HPC SYSTEMS, Inc. All rights 5

87 2 体 :O(2) の計算量 2009/3/ HPC SYSTEMS, Inc. All rights 6

88 逐次版 N 体問題のアルゴリズム m m F= G r b 各質点について (N 回 ) 全ての質点について (N-1 回 ) a c 座標と質量から相互作用の力を計算する d 働く力の総和を計算する次の時刻の座標として座標を更新する 2009/3/ HPC SYSTEMS, Inc. All rights 7

89 N 体 :O(N 2 ) の計算量 2009/3/ HPC SYSTEMS, Inc. All rights 8

90 N 体 : 逐次版ソースコード概要 (1/2) void calculate_forces(void) { int i, j; coords direction; body **b = barray; double distance, magnitude; } O(N 2 ) for (i = 0; i < N; i++) { for (j = 0; j < N; j++) { if (i == j) continue; distance = sqrt( pow((b[i]->p.x - b[j]->p.x), 2.0) + pow((b[i]->p.y - b[j]->p.y), 2.0) ); magnitude = (G * b[i]->m * b[j]->m) / pow(distance, 2.0); direction.x = b[j]->p.x - b[i]->p.x; direction.y = b[j]->p.y - b[i]->p.y; b[i]->f.x = b[i]->f.x + magnitude * direction.x / distance; b[i]->f.y = b[i]->f.y + magnitude * direction.y / distance; } } 反作用の力も同時に求められるがそれを並列化すると隣のプロセスを触ることになり複雑化するため今回は考えない質点間の距離相互作用の力力の総和 2009/3/ HPC SYSTEMS, Inc. All rights 9

91 N 体 : 逐次版ソースコード概要 (2/2) void move_bodies(void) { int i; coords deltav; coords deltap; vk 1 body **b = barray; xk 1 for (i = 0; i < N; i++) { deltav.x = b[i]->f.x / b[i]->m * DT; deltav.y = b[i]->f.y / b[i]->m * DT; deltap.x = (b[i]->v.x + deltav.x / 2) * DT; deltap.y = (b[i]->v.y + deltav.y / 2) * DT; b[i]->v.x = b[i]->v.x + deltav.x; b[i]->v.y = b[i]->v.y + deltav.y; b[i]->p.x = b[i]->p.x + deltap.x; b[i]->p.y = b[i]->p.y + deltap.y; /* Reset force vector */ b[i]->f.x = 0.0; b[i]->f.y = 0.0; } } 第三回 deltav + = vk + fx2 / m1 * t = x + v * t + v + k k k / 2* 次の時刻での速度次の時刻での位置 2009/3/ HPC SYSTEMS, Inc. All rights 10 t

92 N 体問題の並列性時刻 t1: 力の計算 1: a b F= G m m r a b 2 ab 計算 1,2 間に依存性がない互いに独立して実行できる = 並列に計算できる 2009/3/26 力の計算 2: c d 2009 HPC SYSTEMS, Inc. All rights F= G m m r c d 2 cd 11

93 N 体問題の並列解法 1. 質点の配列を root から全プロセスへコピー 2. 各プロセスが担当グループの力と座標を計算 3. タイムステップ更新時に配列を root へ回収再コピー時刻 t1 P P P P 各プロセスは自分の色の質点についてのみ計算する rank /3/ HPC SYSTEMS, Inc. All rights 12

94 N 体 : 並列版の動作 rank 0 例 : 各 CPUはN/4 個の質点について計算する計算通信シミュレーションステップ任意回 t 2009/3/ HPC SYSTEMS, Inc. All rights 13

95 MPI_BODY 型の定義 main() { } MPI 初期化 MPI_BODY 型の登録 read_input(); N の値の共有 bodies[] の確保 mystart[], myend[] の計算計測 1 run(steps); 計測 2 2-1の表示メモリ解放 MPI 終了処理 run(int steps) { N, steps の表示 for (step=0; step<steps; step++) { } } 問題 1 Broadcast calculate_forces(); move_bodies(); 問題 2 Gather ビットマップ出力 2009/3/ HPC SYSTEMS, Inc. All rights 14

96 独自型の定義登録 double 7 個からなる MPI_BODY 型の定義登録 Send/Recvの際要素の型として使えるようになる= 書きやすい! typedef struct { double px, py; double vx, vy; double fx, fy; double m; } Body; // この構造体を単位として Send/Recv するため MPI_Datatype MPI_BODY; MPI_Type_contiguous(7, MPI_DOUBLE, &MPI_BODY); MPI_Type_commit(&MPI_BODY); 2009/3/ HPC SYSTEMS, Inc. All rights 15

97 デモアプリの使い方コンパイル mpicc -o nbody nbody_mpi.c error_wrapper.c Disk.c make も使えます入力./nbody < 入力ファイル > < ステップ数 > random5000.data 100 質点 5000 個ランダム配置出力標準出力シミュレーション全体に何秒かかったかビットマップファイル各シミュレーションステップの質点の座標をプロット./STEP????.bmp 余りは末尾のプロセスが担当 2009/3/ HPC SYSTEMS, Inc. All rights 16

98 問題 1 Broadcast を Send/Recv で書こう // // (1) Broadcast bodies from root process to the other processes // if ( rank == root ) { // // TODO: // Send all bodies to non-root processes // } else { // // TODO: // Receive all bodies from the root process // } rank /3/ HPC SYSTEMS, Inc. All rights 17

99 問題 2 Gather を Send/Recv で書こう // // (2) Gather states of the bodies to the root process // if ( rank == root ) { // // TODO: // Receive each bodies from non-root processes // } else { // // TODO: // Send my bodies to the root process // // HINT: // All bodies are stored as bodies[]. rank 0 // Indices of my bodies are between mystart[rank] and myend[rank]. // Use MPI_BODY as element type. // } 2009/3/ HPC SYSTEMS, Inc. All rights

100 実行時間書き込みページ nproc 実行時間 [ 秒 ] 速度向上率 /3/ HPC SYSTEMS, Inc. All rights 19

101 並列版の性能スケーラビリティ Linear 速度向上率多項式近似速度向上率 /3/ nproc 2009 HPC SYSTEMS, Inc. All rights 20

102 答え合わせ 1 if ( rank == root ) { for ( i = 1; i < nproc; i++ ) { MPI_Send(bodies, N, MPI_BODY, i, 1, MPI_COMM_WORLD); } } else { MPI_Recv(bodies, N, MPI_BODY, root, 1, MPI_COMM_WORLD, &status); } rankが1 以上のプロセスへ Send rankがrootのプロセスから Recv 2009/3/ HPC SYSTEMS, Inc. All rights 21

103 答え合わせ 2 if ( rank == root ) { rankが1 以上の for ( i = 1; i < nproc; i++ ) { プロセスから MPI_Recv(&bodies[myStart[i]], Recv myend[i] - mystart[i], MPI_BODY, i, 2, MPI_COMM_WORLD, &status); } } else { rankがrootの MPI_Send(&bodies[myStart[rank]], プロセスへ myend[rank] - mystart[rank], Send MPI_BODY, root, 2, MPI_COMM_WORLD); } 2009/3/ HPC SYSTEMS, Inc. All rights 22

104 集合通信関数を使うと MPI_Gather(&bodies[myStart[rank]], myend[rank]-mystart[rank], MPI_BODY, &bodies[root], myend[rank]-mystart[rank], root, MPI_COMM_WORLD); rank localc localc localc localc MPI_Send MPI_Send MPI_Send MPI_Send MPI_Recv x4 C 2009/3/ HPC SYSTEMS, Inc. All rights 23

105 MPI 総括並列プログラム開発の標準的 API 基本データを Send/Recv する rank=1 並列性の把握データ分配粒度がカギとにかくネットワーク通信を控えよう 2009/3/26 rank= HPC SYSTEMS, Inc. All rights 1

106 HPCシステムズが提供していくものセミナーでの情報提供メニーコア特殊アーキテクチャ適材適所 Cell, GPU, Vectorized (Specialized) hardware 高速インターコネクト並列向け言語とコンパイラプロファイリング/性能解析並列デバッグメニーコアエキスパート育成 2009/3/26 第三回MPI実践セミナー実践セミナー Copyright (C) 第三回 2009 HPC SYSTEMS, Inc. All rights 2

107 並列化無償代行サービス並列化はやっぱり難しいな誰かやってくれないかな 2009/3/ HPC SYSTEMS, Inc. All rights 3

108 御清聴ありがとうございました質問がありましたらどうぞアンケートご記入をお願いいたしますお帰りの際名札をご返却ください 2009/3/ HPC SYSTEMS, Inc. All rights 4

NUMAの構成

NUMAの構成メッセージパッシングプログラミング天野共有メモリ対メッセージパッシング共有メモリモデル共有変数を用いた単純な記述自動並列化コンパイラ簡単なディレクティブによる並列化 :OpenMP メッセージパッシング形式検証が可能 ( ブロッキング ) 副作用がない ( 共有変数は副作用そのもの ) コストが小さいメッセージパッシングモデル共有変数は使わない共有メモリがないマシンでも実装可能クラスタ

第三回 MPI 実践セミナー HPC システムズ株式会社 新規事業企画室南本和秀 Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 1 rights reserved.

第三回 MPI 実践セミナー HPC システムズ株式会社新規事業企画室南本和秀 Copyright (C) 2009 HPC SYSTEMS, Inc. All 2009/3/26 1 rights reserved.