図4-4-1 領域分割
MPIを使用することにより、計算時間の大部分を占める電磁界の更新を高速化することができます。
ここでは計算領域をX方向にプロセスの数で分割することとします。
図4-4-1は計算領域をX方向に4分割したものです。
ただし、X方向のセル数Nxがプロセス数Npの倍数でないときはなるべく不均一性が小さくなるように分割します。
各プロセスの必要メモリーはプロセス数分の1になります。
これにより計算ノードの数を増やせばそれに比例して計算できる問題サイズの上限が大きくなります。
MPIではプロセス間の通信以外に、
iループのインデックスが0<=i<=NxからiMin<=i<=iMax
と変わるためにプログラムの変更が多数必要になります。
このためにMPIによる高速化はSIMDやマルチスレッドによる高速化に比べて作業量が増えます。
MPIではタイムステップごとに計算領域の境界で隣のプロセスとの通信が必要になります。
通信する成分は図4-4-2の青色すなわちHy,Hz成分です。
プロセスPのHy,HzをプロセスP+1に送信し、
プロセスP+1のHy,HzをプロセスPに送信します(正確には吸収境界条件の電磁界も通信することが必要です)。
通信量は8NyNzバイトと評価されます。
境界線上にある赤色すなわちEy,Ez,Hx成分は両方のプロセスで重複して持ち、
両方のプロセスで更新されます。
リスト4-4-1に電磁界を通信する関数を示します。
通信する配列のアドレス(Offset)と長さを予め計算しバッファを確保しておきます。
変数Offsetの要素0/1/2/3は順に「-X境界の外側/-X境界の内側/+X境界の内側/+X境界の外側」を表します。
プロセスPの-X境界でプロセスP-1とHy,Hzを送受信し、
その後+X境界でプロセスP+1とHy,Hzを送受信します。
通信にはデッドロックの心配のないMPI_Sendrecv関数を使用します。
1 void comm_boundary()
2 {
3 #ifdef MPI
4 MPI_Status status;
5 const int tag = 0;
6 size_t size_hy = Length_Hy * sizeof(float);
7 size_t size_hz = Length_Hz * sizeof(float);
8 int count = Length_Hy + Length_Hz;
9
10 // -X boundary
11 if (commRank > 0) {
12 // copy to buffer
13 memcpy(sendBuf, Hy + Offset_Hy[1], size_hy);
14 memcpy(sendBuf + Length_Hy, Hz + Offset_Hz[1], size_hz);
15
16 // MPI
17 int dst = commRank - 1;
18 MPI_Sendrecv(sendBuf, count, MPI_FLOAT, dst, tag,
19 recvBuf, count, MPI_FLOAT, dst, tag, MPI_COMM_WORLD, &status);
20
21 // copy from buffer
22 memcpy(Hy + Offset_Hy[0], recvBuf, size_hy);
23 memcpy(Hz + Offset_Hz[0], recvBuf + Length_Hy, size_hz);
24 }
25
26 // +X boundary
27 if (commRank < commSize - 1) {
28 // copy to buffer
29 memcpy(sendBuf, Hy + Offset_Hy[2], size_hy);
30 memcpy(sendBuf + Length_Hy, Hz + Offset_Hz[2], size_hz);
31
32 // MPI
33 int dst = commRank + 1;
34 MPI_Sendrecv(sendBuf, count, MPI_FLOAT, dst, tag,
35 recvBuf, count, MPI_FLOAT, dst, tag, MPI_COMM_WORLD, &status);
36
37 // copy from buffer
38 memcpy(Hy + Offset_Hy[3], recvBuf, size_hy);
39 memcpy(Hz + Offset_Hz[3], recvBuf + Length_Hy, size_hz);
40 }
41 #endif
42 }
収束判定に用いる平均電磁界については各プロセスの部分和の和を計算し、
結果を全プロセスで共有し全プロセスで収束判定を行います。
リスト4-4-2にソースコードを示します。
上記の作業を行うためにMPI_Allreduce関数を使用しています。
1 void comm_average(double fsum[])
2 {
3 #ifdef MPI
4 double ftmp[2];
5
6 MPI_Allreduce(fsum, ftmp, 2, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
7
8 fsum[0] = ftmp[0];
9 fsum[1] = ftmp[1];
10 #endif
11 }
アンテナの入力インピーダンスを求めるには、反復計算が終了した後、
給電点の電圧と電流の時間波形をルートプロセスに集めることが必要です。
リスト4-4-3にソースコードを示します。
給電点が非ルートプロセスの領域にあるときだけ通信が必要になります。
給電点を有しているプロセスが時間波形をルートプロセスに送信し、
ルートプロセスがそれを受信します。
1 void comm_feed()
2 {
3 #ifdef MPI
4 MPI_Status status;
5
6 int count = Solver.maxiter + 1;
7 for (int n = 0; n < NFeed; n++) {
8 // non-root only
9 if ((commRank == 0) && !iProc(Feed[n].i)) {
10 MPI_Recv(&VFeed[n * count], count, MPI_FLOAT, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
11 MPI_Recv(&IFeed[n * count], count, MPI_FLOAT, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, &status);
12 }
13 else if ((commRank > 0) && iProc(Feed[n].i)) {
14 MPI_Send(&VFeed[n * count], count, MPI_FLOAT, 0, 0, MPI_COMM_WORLD);
15 MPI_Send(&IFeed[n * count], count, MPI_FLOAT, 0, 0, MPI_COMM_WORLD);
16 }
17 MPI_Barrier(MPI_COMM_WORLD);
18 }
19 #endif
20 }
ここでは近傍界のうち観測面の通信について述べます。
観測面にはX/Y/Z面の3通りがあり、
かつ領域がX方向に分割されているために処理が煩雑になります。
反復計算を行う前に準備作業として各プロセスのデータ数を計算しておきます。
リスト4-4-4がそのための関数です。
ソースコードの説明
5行目:各プロセスの近傍界観測面ごとのデータの数を保存する配列を作成します。
7-40行目:その配列に値を代入します。
14,24,34行目:インデックスiを有するプロセスのみカウントします。
"iProc(i)"はインデックスiが属するプロセス番号です。
1 static void setupSizeNear2d()
2 {
3 if (NNear2d <= 0) return;
4
5 l_LNear2d = (int *)malloc(NNear2d * sizeof(int));
6
7 for (int m = 0; m < NNear2d; m++) {
8 int num = 0;
9 if (Near2d[m].dir == 'X') {
10 // Y-Z
11 int i = Near2d[m].id0;
12 for (int j = 0; j < Ny; j++) {
13 for (int k = 0; k < Nz; k++) {
14 if (!iProc(i)) continue;
15 num++;
16 }
17 }
18 }
19 else if (Near2d[m].dir == 'Y') {
20 // X-Z
21 //int j = Near2d[m].id0;
22 for (int i = 0; i < Nx; i++) {
23 for (int k = 0; k < Nz; k++) {
24 if (!iProc(i)) continue;
25 num++;
26 }
27 }
28 }
29 else if (Near2d[m].dir == 'Z') {
30 // X-Y
31 //int k = Near2d[m].id0;
32 for (int i = 0; i < Nx; i++) {
33 for (int j = 0; j < Ny; j++) {
34 if (!iProc(i)) continue;
35 num++;
36 }
37 }
38 }
39 l_LNear2d[m] = num;
40 }
41
42 int sum = 0;
43 for (int m = 0; m < NNear2d; m++) {
44 sum += l_LNear2d[m];
45 }
46 Near2d_size = sum * NFreq2 * sizeof(f_complex_t);
47 }
反復計算が終了した後に、
各プロセスの持っている近傍界をルートプロセスに集めることが必要です。
リスト4-4-5がそのための関数です。
通信するデータの大きさがプロセスによって異なりますので通信にはMPI_Gatherv関数を使用します。
ソースコードの説明
11-31行目:ルートプロセスが観測面全体の配列を確保します。
34-37行目:ルートプロセスが各プロセスの持っているデータが全体配列の中でどの位置にあるかを表すoffset,count配列を確保します。
50-51行目:各プロセスの持っているデータの大きさをルートプロセスのcount配列に集めます。
54-59行目:ルートプロセスがcount配列を元にoffset配列を計算します。
62-63行目:全プロセスが送信用バッファを用意します。
66-69行目:ルートプロセスが受信用バッファを用意します。
86-124行目:ルートプロセスに近傍界を集めます。
MPI_Gatherv関数の引数に上で計算したcount,offsetが含まれています。1周波数につき1回通信しています。
1 void comm_near2d()
2 {
3 #ifdef MPI
4 int id;
5 int g_num = 0;
6 int64_t g_adr = 0;
7 int *offset = NULL, *count = NULL;
8 float *send = NULL, *recv = NULL;
9
10 // alloc global array (root)
11 if (commRank == 0) {
12 int sum = 0;
13 for (int m = 0; m < NNear2d; m++) {
14 if (Near2d[m].dir == 'X') {
15 sum += Ny * Nz;
16 }
17 else if (Near2d[m].dir == 'Y') {
18 sum += Nz * Nx;
19 }
20 else if (Near2d[m].dir == 'Z') {
21 sum += Nx * Ny;
22 }
23 }
24 size_t g_size = sum * NFreq2 * sizeof(f_complex_t);
25 Near2dEx = (f_complex_t *)malloc(g_size);
26 Near2dEy = (f_complex_t *)malloc(g_size);
27 Near2dEz = (f_complex_t *)malloc(g_size);
28 Near2dHx = (f_complex_t *)malloc(g_size);
29 Near2dHy = (f_complex_t *)malloc(g_size);
30 Near2dHz = (f_complex_t *)malloc(g_size);
31 }
32
33 // alloc offset and count (root)
34 if (commRank == 0) {
35 offset = (int *)malloc(commSize * sizeof(int));
36 count = (int *)malloc(commSize * sizeof(int));
37 }
38
39 int64_t l_adr = 0;
40 if (commRank == 0) {
41 g_adr = 0;
42 }
43 for (int m = 0; m < NNear2d; m++) {
44 // local number of data
45 int l_num = l_LNear2d[m];
46
47 // gather count to root
48 // 6 = Ex/Ey/Ez/Hx/Hy/Hz components
49 // 2 = sizeof(f_complex_t) / sizeof(float)
50 int l_count = l_num * 6 * 2;
51 MPI_Gather(&l_count, 1, MPI_INT, count, 1, MPI_INT, 0, MPI_COMM_WORLD);
52
53 // setup offset and count (root)
54 if (commRank == 0) {
55 offset[0] = 0;
56 for (int iproc = 1; iproc < commSize; iproc++) {
57 offset[iproc] = offset[iproc - 1] + count[iproc - 1];
58 }
59 }
60
61 // send buffer (all process)
62 size_t size = l_count * sizeof(float);
63 send = (float *)malloc(size);
64
65 // recv buffer (global array : root only)
66 if (commRank == 0) {
67 size = (offset[commSize - 1] + count[commSize - 1]) * sizeof(float);
68 recv = (float *)malloc(size);
69 }
70
71 // global number of data
72 if (commRank == 0) {
73 g_num = 0;
74 if (Near2d[m].dir == 'X') {
75 g_num = Ny * Nz;
76 }
77 else if (Near2d[m].dir == 'Y') {
78 g_num = Nz * Nx;
79 }
80 else if (Near2d[m].dir == 'Z') {
81 g_num = Nx * Ny;
82 }
83 }
84
85 // gather E and H to root
86 for (int ifreq = 0; ifreq < NFreq2; ifreq++) {
87 id = 0;
88 for (int n = 0; n < l_num; n++) {
89 send[id++] = l_Near2dEx[l_adr].r;
90 send[id++] = l_Near2dEx[l_adr].i;
91 send[id++] = l_Near2dEy[l_adr].r;
92 send[id++] = l_Near2dEy[l_adr].i;
93 send[id++] = l_Near2dEz[l_adr].r;
94 send[id++] = l_Near2dEz[l_adr].i;
95 send[id++] = l_Near2dHx[l_adr].r;
96 send[id++] = l_Near2dHx[l_adr].i;
97 send[id++] = l_Near2dHy[l_adr].r;
98 send[id++] = l_Near2dHy[l_adr].i;
99 send[id++] = l_Near2dHz[l_adr].r;
100 send[id++] = l_Near2dHz[l_adr].i;
101 l_adr++;
102 }
103
104 MPI_Gatherv(send, id, MPI_FLOAT, recv, count, offset, MPI_FLOAT, 0, MPI_COMM_WORLD);
105
106 if (commRank == 0) {
107 id = 0;
108 for (int n = 0; n < g_num; n++) {
109 Near2dEx[g_adr].r = recv[id++];
110 Near2dEx[g_adr].i = recv[id++];
111 Near2dEy[g_adr].r = recv[id++];
112 Near2dEy[g_adr].i = recv[id++];
113 Near2dEz[g_adr].r = recv[id++];
114 Near2dEz[g_adr].i = recv[id++];
115 Near2dHx[g_adr].r = recv[id++];
116 Near2dHx[g_adr].i = recv[id++];
117 Near2dHy[g_adr].r = recv[id++];
118 Near2dHy[g_adr].i = recv[id++];
119 Near2dHz[g_adr].r = recv[id++];
120 Near2dHz[g_adr].i = recv[id++];
121 g_adr++;
122 }
123 }
124 }
125
126 // free
127 free(send);
128 if (commRank == 0) {
129 free(recv);
130 }
131 }
132
133 // free
134 if (commRank == 0) {
135 free(offset);
136 free(count);
137 }
138 #endif
139 }
表4-4-1にWindows環境でプロセス数を1〜8と変えたときの計算時間の測定結果を示します。
2プロセスである程度速くなりますが、プロセス数をそれ以上増やしても速くなりません。
また2スレッド以上ではSIMDの有無で計算時間が変わりません。
これはOpenMPと同じ結果です。
| プロセス数 | ベンチマークNo.1 | ベンチマークNo.2 | ||
|---|---|---|---|---|
| SIMDなし | SSE | SIMDなし | SSE | |
| 1 | 28.9秒(1.0) | 20.0秒(1.0) | 233.6秒(1.0) | 156.7秒(1.0) |
| 2 | 18.0秒(1.61) | 17.0秒(1.18) | 141.3秒(1.65) | 129.8秒(1.21) |
| 4 | 17.9秒(1.61) | 17.8秒(1.12) | 134.4秒(1.74) | 129.4秒(1.21) |
| 8 | 18.8秒(1.54) | 18.7秒(1.07) | 137.8秒(1.70) | 135.6秒(1.16) |
表4-4-2と図4-4-3にLinux環境でプロセス数を1〜40と変えたときの計算時間の測定結果を示します。
20プロセスまで速くなり、これはOpenMPより良好な結果です。
| プロセス数 | ベンチマークNo.1 | ベンチマークNo.2 | ||
|---|---|---|---|---|
| SIMDなし | SSE | SIMDなし | SSE | |
| 1 | 31.4秒(1.0) | 19.2秒(1.0) | 251.9秒(1.0) | 188.1秒(1.0) |
| 2 | 14.9秒(2.11) | 8.1秒(2.37) | 115.9秒(2.17) | 69.8秒(2.69) |
| 4 | 7.9秒(3.97) | 4.4秒(4.36) | 63.2秒(3.99) | 43.9秒(4.28) |
| 10 | 3.9秒(8.05) | 3.0秒(6.40) | 42.7秒(5.90) | 41.4秒(4.54) |
| 20 | 2.3秒(13.65) | 1.7秒(11.29) | 28.1秒(8.96) | 28.1秒(6.69) |
| 40 | 6.3秒(4.98) | 4.0秒(4.80) | 35.2秒(7.16) | 35.7秒(5.27) |
MPIは複数台のコンピュータで並列計算することができ、これが本来のMPIの目的です。
複数台で計算するための方法については[3]を参考にしてください。
MPIでは1台のコンピュータをノードと呼びます。
表4-4-3に1〜5ノードの計算時間を示します。
CPUは8コア16スレッドであり、8コアを使用しています。
ほぼノード数に比例して速くなることがわかります。
前節のようにコア数の多い1ノードで計算するよりも、
本節のように各ノードのコア数が少なくても複数のノードで計算したほうが速くなることがあります。
| ノード数 | プロセス数 | ベンチマークNo.2(SSE) |
|---|---|---|
| 1 | 8 | 66.1秒(1.0) |
| 2 | 16 | 36.1秒(1.83) |
| 3 | 24 | 26.6秒(2.48) |
| 4 | 32 | 21.8秒(3.03) |
| 5 | 40 | 17.1秒(3.87) |