4.4 MPIによる高速化

4.4.1 領域分割

MPIを使用することにより、計算時間の大部分を占める電磁界の更新を高速化することができます。

ここでは計算領域をX方向にプロセスの数で分割することとします。
図4-4-1は計算領域をX方向に4分割したものです。 ただし、X方向のセル数Nxがプロセス数Npの倍数でないときはなるべく不均一性が小さくなるように分割します。
各プロセスの必要メモリーはプロセス数分の1になります。 これにより計算ノードの数を増やせばそれに比例して計算できる問題サイズの上限が大きくなります。

MPIではプロセス間の通信以外に、 iループのインデックスが0<=i<=NxからiMin<=i<=iMax と変わるためにプログラムの変更が多数必要になります。
このためにMPIによる高速化はSIMDやマルチスレッドによる高速化に比べて作業量が増えます。


図4-4-1 領域分割

4.4.2 電磁界の通信

MPIではタイムステップごとに計算領域の境界で隣のプロセスとの通信が必要になります。
通信する成分は図4-4-2の青色すなわちHy,Hz成分です。 プロセスPのHy,HzをプロセスP+1に送信し、 プロセスP+1のHy,HzをプロセスPに送信します(正確には吸収境界条件の電磁界も通信することが必要です)。
通信量は8NyNzバイトと評価されます。
境界線上にある赤色すなわちEy,Ez,Hx成分は両方のプロセスで重複して持ち、 両方のプロセスで更新されます。


図4-4-2 電磁界の通信

リスト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関数を使用します。


リスト4-4-1 電磁界の通信
     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.3 平均電磁界の通信

収束判定に用いる平均電磁界については各プロセスの部分和の和を計算し、 結果を全プロセスで共有し全プロセスで収束判定を行います。
リスト4-4-2にソースコードを示します。 上記の作業を行うためにMPI_Allreduce関数を使用しています。


リスト4-4-2 平均電磁界の通信
     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.4 時間波形の通信

アンテナの入力インピーダンスを求めるには、反復計算が終了した後、 給電点の電圧と電流の時間波形をルートプロセスに集めることが必要です。
リスト4-4-3にソースコードを示します。
給電点が非ルートプロセスの領域にあるときだけ通信が必要になります。 給電点を有しているプロセスが時間波形をルートプロセスに送信し、 ルートプロセスがそれを受信します。


リスト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	}

4.4.5 近傍界の通信

ここでは近傍界のうち観測面の通信について述べます。
観測面にはX/Y/Z面の3通りがあり、 かつ領域がX方向に分割されているために処理が煩雑になります。
反復計算を行う前に準備作業として各プロセスのデータ数を計算しておきます。
リスト4-4-4がそのための関数です。

ソースコードの説明
5行目:各プロセスの近傍界観測面ごとのデータの数を保存する配列を作成します。
7-40行目:その配列に値を代入します。
14,24,34行目:インデックスiを有するプロセスのみカウントします。 "iProc(i)"はインデックスiが属するプロセス番号です。


リスト4-4-4 近傍界の管理
     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回通信しています。


リスト4-4-5 近傍界の通信
     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.6 MPIの計算時間

表4-4-1にWindows環境でプロセス数を1〜8と変えたときの計算時間の測定結果を示します。
2プロセスである程度速くなりますが、プロセス数をそれ以上増やしても速くなりません。
また2スレッド以上ではSIMDの有無で計算時間が変わりません。
これはOpenMPと同じ結果です。

表4-4-1 MPIの計算時間(Windows環境、i7-4770K、()内は1プロセスとの速度比)
プロセス数ベンチマークNo.1ベンチマークNo.2
SIMDなしSSESIMDなしSSE
128.9秒(1.0)20.0秒(1.0)233.6秒(1.0)156.7秒(1.0)
218.0秒(1.61)17.0秒(1.18)141.3秒(1.65)129.8秒(1.21)
417.9秒(1.61)17.8秒(1.12)134.4秒(1.74)129.4秒(1.21)
818.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より良好な結果です。

表4-4-2 MPIの計算時間(FOCUSスパコンDシステム、Xeon E5-2670v2 x2、icc17.0、()内は1プロセスとの速度比)
プロセス数ベンチマークNo.1ベンチマークNo.2
SIMDなしSSESIMDなしSSE
131.4秒(1.0)19.2秒(1.0)251.9秒(1.0)188.1秒(1.0)
214.9秒(2.11)8.1秒(2.37)115.9秒(2.17)69.8秒(2.69)
47.9秒(3.97)4.4秒(4.36)63.2秒(3.99)43.9秒(4.28)
103.9秒(8.05)3.0秒(6.40)42.7秒(5.90)41.4秒(4.54)
202.3秒(13.65)1.7秒(11.29)28.1秒(8.96)28.1秒(6.69)
406.3秒(4.98)4.0秒(4.80)35.2秒(7.16)35.7秒(5.27)


図4-4-3 MPIの計算時間(FOCUSスパコンDシステム、Xeon E5-2670v2 x2、icc17.0)

4.4.7 複数台での並列計算

MPIは複数台のコンピュータで並列計算することができ、これが本来のMPIの目的です。
複数台で計算するための方法については[3]を参考にしてください。
MPIでは1台のコンピュータをノードと呼びます。
表4-4-3に1〜5ノードの計算時間を示します。 CPUは8コア16スレッドであり、8コアを使用しています。
ほぼノード数に比例して速くなることがわかります。
前節のようにコア数の多い1ノードで計算するよりも、 本節のように各ノードのコア数が少なくても複数のノードで計算したほうが速くなることがあります。

表4-4-3 複数台の計算時間(FOCUSスパコンHシステム、Xeon D-1541 x1、icc17.0、10Gbps、()内は1ノードとの速度比)
ノード数プロセス数ベンチマークNo.2(SSE)
1866.1秒(1.0)
21636.1秒(1.83)
32426.6秒(2.48)
43221.8秒(3.03)
54017.1秒(3.87)