2.4 並列計算

2.4.1 OpenMPを用いたCPUの並列計算

リスト2-4-1にCPU版のred-black SOR法の反復計算部分を示します。
13行目でOpenMPを用いてループiを並列化しています。 またreduction句を用いて残差の和を求めています。[6]
19行目でred-black法の判定を行っています。
電極は配列"idVolt"によって間接的に参照しています。idVolt=0は非電極(すなわち誘電体)、 idVolt=1,2,...は電圧番号1,2,...を表します。 idVolt>0のときは電圧は既知であるために計算を飛ばしています。
誘電体は配列"idEpsr"によって間接的に参照しています。idEpsr=0は空気(比誘電率=1)、 idEpsr=1,2,...は比誘電率番号1,2,...を表します。

リスト2-4-1 OpenMPを用いた反復計算部の並列化(CPU版)
     1	static real_t update_rb(int oe)
     2	{
     3		real_t res2 = 0;
     4	
     5		real_t omega = (real_t)Solver.omega;
     6	
     7		int64_t ni = D(1, 0, 0) - D(0, 0, 0);
     8		int64_t nj = D(0, 1, 0) - D(0, 0, 0);
     9		int64_t nk = D(0, 0, 1) - D(0, 0, 0);
    10	
    11		int i;
    12	#ifdef _OPENMP
    13	#pragma omp parallel for reduction (+: res2)
    14	#endif
    15		for (    i = 0; i <= Nx; i++) {
    16		for (int j = 0; j <= Ny; j++) {
    17		for (int k = 0; k <= Nz; k++) {
    18			int64_t n = D(i, j, k);
    19			if (((i + j + k) % 2 == oe) && !idVolt[n]) {
    20				if      (i == 0) {
    21					V[n] = V[n + ni];
    22				}
    23				else if (i == Nx) {
    24					V[n] = V[n - ni];
    25				}
    26				else if (j == 0) {
    27					V[n] = V[n + nj];
    28				}
    29				else if (j == Ny) {
    30					V[n] = V[n - nj];
    31				}
    32				else if (k == 0) {
    33					V[n] = V[n + nk];
    34				}
    35				else if (k == Nz) {
    36					V[n] = V[n - nk];
    37				}
    38				else {
    39					real_t e = Epsr[idEpsr[n]];
    40					real_t axp = (e + Epsr[idEpsr[n + ni]]) * RXp[i];
    41					real_t axm = (e + Epsr[idEpsr[n - ni]]) * RXm[i];
    42					real_t ayp = (e + Epsr[idEpsr[n + nj]]) * RYp[j];
    43					real_t aym = (e + Epsr[idEpsr[n - nj]]) * RYm[j];
    44					real_t azp = (e + Epsr[idEpsr[n + nk]]) * RZp[k];
    45					real_t azm = (e + Epsr[idEpsr[n - nk]]) * RZm[k];
    46					real_t res = omega * ((
    47						axp * V[n + ni] +
    48						axm * V[n - ni] +
    49						ayp * V[n + nj] +
    50						aym * V[n - nj] +
    51						azp * V[n + nk] +
    52						azm * V[n - nk]) / (axp + axm + ayp + aym + azp + azm) - V[n]);
    53					V[n] += res;
    54					res2 += res * res;
    55				}
    56			}
    57		}
    58		}
    59		}
    60	
    61		return res2;
    62	}

2.4.2 CUDAを用いたGPUの並列計算

リスト2-4-2にGPU版のred-black SOR法の反復計算部分を示します。
blockは(k, j)の2次元の(16, 16)とし、gridは((Nz+1)/16, (Ny+1)/16, Nx+1)の3次元とします。
これは(i, j, k)の3次元空間を1次元配列で並べるときに内側からk, j, iの順にとっているためです。 CUDAでは変数のアドレスがスレッド番号の順に並んでいることが重要です。[8]
残差の計算はblock内で部分和をshared memoryのreduction演算で求め、 部分和の配列をCPUに転送しその和をCPUで計算しています。 このようにすると変数全体を転送して残差を計算することに比べると転送量が1/256倍になり、 データ転送に要する時間を無視することができます。

リスト2-4-2 CUDAを用いた反復計算部の並列化(GPU版)
     1	__host__ __device__ real_t update(
     2		int i, int j, int k, int nx, int ny, int nz, int64_t ni, int64_t nj, int64_t nk, int64_t n,
     3		real_t *rxp, real_t *ryp, real_t *rzp, real_t *rxm, real_t *rym, real_t *rzm,
     4		real_t *v, unsigned char *idepsr, real_t *epsr, real_t omega)
     5	{
     6		real_t res = 0;
     7	
     8		if      (i == 0) {
     9			v[n] = v[n + ni];
    10		}
    11		else if (i == nx) {
    12			v[n] = v[n - ni];
    13		}
    14		else if (j == 0) {
    15			v[n] = v[n + nj];
    16		}
    17		else if (j == ny) {
    18			v[n] = v[n - nj];
    19		}
    20		else if (k == 0) {
    21			v[n] = v[n + nk];
    22		}
    23		else if (k == nz) {
    24			v[n] = v[n - nk];
    25		}
    26		else {
    27			real_t e = epsr[idepsr[n]];
    28	
    29			real_t axp = (e + epsr[idepsr[n + ni]]) * rxp[i];
    30			real_t axm = (e + epsr[idepsr[n - ni]]) * rxm[i];
    31			real_t ayp = (e + epsr[idepsr[n + nj]]) * ryp[j];
    32			real_t aym = (e + epsr[idepsr[n - nj]]) * rym[j];
    33			real_t azp = (e + epsr[idepsr[n + nk]]) * rzp[k];
    34			real_t azm = (e + epsr[idepsr[n - nk]]) * rzm[k];
    35	
    36			res = omega * ((
    37				axp * v[n + ni] +
    38				axm * v[n - ni] +
    39				ayp * v[n + nj] +
    40				aym * v[n - nj] +
    41				azp * v[n + nk] +
    42				azm * v[n - nk]) / (axp + axm + ayp + aym + azp + azm) - v[n]);
    43			v[n] += res;
    44		}
    45	
    46		return (res * res);
    47	}
    48	
    49	// GPU
    50	__global__
    51	static void update_gpu(int oe,
    52		int nx, int ny, int nz, int64_t ni, int64_t nj, int64_t nk,
    53		real_t *rxp, real_t *ryp, real_t *rzp, real_t *rxm, real_t *rym, real_t *rzm,
    54		real_t *v, unsigned char *idvolt, unsigned char *idepsr, real_t *epsr, real_t omega, real_t *res2)
    55	{
    56		extern __shared__ real_t sm[];
    57	
    58		int i = blockIdx.z;
    59		int j = threadIdx.y + (blockIdx.y * blockDim.y);
    60		int k = threadIdx.x + (blockIdx.x * blockDim.x);
    61	
    62		int nthread = blockDim.x * blockDim.y;
    63		int tid = threadIdx.x + (threadIdx.y * blockDim.x);
    64		int bid = blockIdx.x + (blockIdx.y * gridDim.x) + (blockIdx.z * gridDim.x * gridDim.y);
    65	
    66		int64_t n = (i * ni) + (j * nj) + (k * nk);
    67	
    68		if ((i <= nx) &&
    69		    (j <= ny) &&
    70		    (k <= nz) &&
    71		    ((i + j + k) % 2 == oe) && !idvolt[n]) {
    72			sm[tid] = update(
    73				i, j, k, nx, ny, nz, ni, nj, nk, n,
    74				rxp, ryp, rzp, rxm, rym, rzm,
    75				v, idepsr, epsr, omega);
    76		}
    77		else {
    78			sm[tid] = 0;
    79		}
    80	
    81		reduction_sum(tid, nthread, sm, &res2[bid]);
    82	}

2.4.3 MPIを用いた並列計算

MPIを用いて並列計算するには図2-4-1のように計算領域を分割します。 本シミュレーターではX方向に等分割します。
プロセス数をPとすると計算時間が1/Pになります。
MPIによる並列計算はCPU版とGPU版の両方に適用することができます。
各プロセスの必要メモリーについては3.2で述べるように全領域の電圧分布をファイルに保存するときは逐次版と同じですが、 ファイルに保存しないときは1/Pになります。 ただしGPU版ではファイルの保存はCPUで行いますので、GPUの必要メモリーは常に1/Pです。


図2-4-1 計算領域のプロセスへの分割(4プロセスの場合)

2.1で述べたように差分法では隣の節点の電圧値を使用しますので、 各プロセスのX方向の両側の境界面では隣のプロセスの電圧値が必要になります。 従ってMPIの通信によって送信と受信を行います。 通信は各反復計算ごとに必要になります。red-black法では反復計算ごとに2回必要になります。 通信量は(Ny+1)*(Nz+1)ワードですのでNx/Pが十分大きければ通信に要する時間は無視することができます。


図2-4-2 隣接プロセスとの通信

リスト2-4-3に通信部を示します。 送受信時にデッドロックが発生しないように設計されているMPI_Sendrecv関数を使用します。 -X方向の境界面と+X方向の境界面の2回必要になります。

リスト2-4-3 隣接プロセスとの送受信
     1	void comm_boundary(real_t *sendbuf, real_t *recvbuf)
     2	{
     3	#ifdef MPI
     4		MPI_Status status;
     5		const int tag = 0;
     6		const int count = (Ny + 1) * (Nz + 1);
     7	
     8		// -X boundary
     9		if (commRank > 0) {
    10			// copy to buffer
    11			int isend = 0;
    12			for (int j = 0; j <= Ny; j++) {
    13			for (int k = 0; k <= Nz; k++) {
    14				sendbuf[isend++] = V[D(1, j, k)];
    15			}
    16			}
    17			assert(isend == count);
    18	
    19			// MPI
    20			int dst = commRank - 1;
    21			MPI_Sendrecv(sendbuf, count, MPI_REAL_T, dst, tag,
    22			             recvbuf, count, MPI_REAL_T, dst, tag, MPI_COMM_WORLD, &status);
    23	
    24			// copy from buffer
    25			int irecv = 0;
    26			for (int j = 0; j <= Ny; j++) {
    27			for (int k = 0; k <= Nz; k++) {
    28				V[D(0, j, k)] = recvbuf[irecv++];
    29			}
    30			}
    31			assert(irecv == count);
    32		}
    33	
    34		// +X boundary
    35		if (commRank < commSize - 1) {
    36			// copy to buffer
    37			int isend = 0;
    38			for (int j = 0; j <= Ny; j++) {
    39			for (int k = 0; k <= Nz; k++) {
    40				sendbuf[isend++] = V[D(iMax2 - iMin2 - 1, j, k)];
    41			}
    42			}
    43			assert(isend == count);
    44	
    45			// MPI
    46			int dst = commRank + 1;
    47			MPI_Sendrecv(sendbuf, count, MPI_REAL_T, dst, tag,
    48			             recvbuf, count, MPI_REAL_T, dst, tag, MPI_COMM_WORLD, &status);
    49	
    50			// copy from buffer
    51			int irecv = 0;
    52			for (int j = 0; j <= Ny; j++) {
    53			for (int k = 0; k <= Nz; k++) {
    54				V[D(iMax2 - iMin2, j, k)] = recvbuf[irecv++];
    55			}
    56			}
    57			assert(irecv == count);
    58		}
    59	#endif
    60	}