3.4 MPIによる高速化

3.4.1 領域分割

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

図3-4-1 領域分割

3.4.2 電界の通信

MPIでは行列・ベクトル積の計算の前に計算領域の境界で隣のプロセスとの通信が必要になります。
プロセスPは領域境界より左の成分を所有し、プロセスP+1は領域境界上とその右の成分を所有します。
通信する成分は図3-4-1に矢印で示した成分です。 このうち、赤い成分はプロセスPからプロセスP+1の方向に通信し、 青い成分はプロセスP+1からプロセスPの方向に通信します。
これによりプロセスP,P+1ともに自分の所有する成分を更新することができます。

図3-4-2 電界の通信

リスト3-4-1にMPIにより電界を通信するプログラムを示します。
アドレスと通信量は予め計算し、バッファーも予め確保しています。

void comm_boundary(d_complex_t *v)
{
#ifdef MPI
	MPI_Status status;
	const int tag = 0;

	// -X direction
	if (commRank > 0) {
		// copy to send buffer
		memcpy(Sendcomplex_m, &v[Send0_m], Sendcount_m * sizeof(d_complex_t));
		for (int64_t n = 0; n < Sendcount_m; n++) {
			Sendbuf_m[(2 * n) + 0] = Sendcomplex_m[n].r;
			Sendbuf_m[(2 * n) + 1] = Sendcomplex_m[n].i;
		}
		// MPI
		int dst = commRank - 1;
		MPI_Sendrecv(Sendbuf_m, 2 * Sendcount_m, MPI_DOUBLE, dst, tag,
		             Recvbuf_m, 2 * Recvcount_m, MPI_DOUBLE, dst, tag, MPI_COMM_WORLD, &status);
		// copy from recv buffer
		for (int64_t n = 0; n < Recvcount_m; n++) {
			 Recvcomplex_m[n].r = Recvbuf_m[(2 * n) + 0];
			 Recvcomplex_m[n].i = Recvbuf_m[(2 * n) + 1];
		}
		memcpy(&v[Recv0_m], Recvcomplex_m, Recvcount_m * sizeof(d_complex_t));
	}

	// +X direction
	if (commRank < commSize - 1) {
		// copy to send buffer
		memcpy(Sendcomplex_p, &v[Send0_p], Sendcount_p * sizeof(d_complex_t));
		for (int64_t n = 0; n < Sendcount_p; n++) {
			Sendbuf_p[(2 * n) + 0] = Sendcomplex_p[n].r;
			Sendbuf_p[(2 * n) + 1] = Sendcomplex_p[n].i;
		}
		// MPI
		int dst = commRank + 1;
		MPI_Sendrecv(Sendbuf_p, 2 * Sendcount_p, MPI_DOUBLE, dst, tag,
		             Recvbuf_p, 2 * Recvcount_p, MPI_DOUBLE, dst, tag, MPI_COMM_WORLD, &status);
		// copy from recv buffer
		for (int64_t n = 0; n < Recvcount_p; n++) {
			 Recvcomplex_p[n].r = Recvbuf_p[(2 * n) + 0];
			 Recvcomplex_p[n].i = Recvbuf_p[(2 * n) + 1];
		}
		memcpy(&v[Recv0_p], Recvcomplex_p, Recvcount_p * sizeof(d_complex_t));
	}
#endif
}

3.4.3 MPIの計算時間

表3-4-1と図3-4-3にプロセス数と計算時間の関係を示します。
OpenMPと同じくプロセス数が少ないときはmatrixモードの方が速いですが、 40プロセスでほぼ同じになります。

図3-4-3 MPIの計算時間