4.4 CUDAによる高速化

4.4.1 CUDAプログラムの設計指針

CUDAによるFDTD法の高速化では下記の点に注意してプログラムを作成します。

  1. 電磁界の更新は並列計算に適したアルゴリズムですので、自然な形で実装できます。
  2. 電磁界の各成分の更新ごとに1回のカーネルを起動します。
  3. blockは(Z,Y)方向の2次元、gridは(Z,Y,X)方向の3次元とします。

4.4.2 電磁界の更新

リスト4-4-1にEx成分を更新する関数を示します。
[3]の手法を用いてCPUコードとGPUコードを共通化しています。
なお、MPIとソースコードを共通化するためにX方向に領域分割した記述にしています。

ソースコードの説明
1行目:計算の実体である関数はCPUとGPUで共通です。 従って関数修飾子"__host__ __device__"を付けます。
12行目:GPUコードには関数修飾子"__global__"を付けます。
18-20行目:gridとblockからi,j,kを計算します。
21-23行目:i,j,kが範囲内にあるときだけ計算を行います(セル数がブロックサイズの倍数とは限らないため)。
53-56行目:blockは内側から順にk,j,iとします。 アドレスはkについて連続ですのでこのようにする必要があります(coalescing of memory)。 CEILは繰り上げ商を計算する自作マクロでありCUDAでは頻繁に使われます。
61行目:unified memoryでは最後に同期を取っています。


リスト4-4-1 Ex成分の更新(updateEx.cu)
     1	__host__ __device__
     2	static void updateEx_(
     3		int64_t n, int nj, int nk,
     4		float *ex, float *hy, float *hz,
     5		float c1, float c2, float ryn, float rzn)
     6	{
     7		ex[n] = c1 * ex[n]
     8		      + c2 * (ryn * (hz[n] - hz[n - nj])
     9		            - rzn * (hy[n] - hy[n - nk]));
    10	}
    11	
    12	__global__
    13	static void updateEx_gpu(
    14		int imin, int jmin, int kmin, int imax, int jmax, int kmax, int ni, int nj, int nk, int n0,
    15		float *ex, float *hy, float *hz, unsigned char *iex,
    16		float *c1, float *c2, float *ryn, float *rzn)
    17	{
    18		int i = imin + threadIdx.z + (blockIdx.z * blockDim.z);
    19		int j = jmin + threadIdx.y + (blockIdx.y * blockDim.y);
    20		int k = kmin + threadIdx.x + (blockIdx.x * blockDim.x);
    21		if ((i <  imax) &&
    22		    (j <= jmax) &&
    23		    (k <= kmax)) {
    24			int64_t n = (i * ni) + (j * nj) + (k * nk) + n0;
    25			updateEx_(
    26				n, nj, nk,
    27				ex, hy, hz,
    28				c1[iex[n]], c2[iex[n]], ryn[j], rzn[k]);
    29		}
    30	}
    31	
    32	static void updateEx_cpu(
    33		int imin, int jmin, int kmin, int imax, int jmax, int kmax, int ni, int nj, int nk, int n0,
    34		float *ex, float *hy, float *hz, unsigned char *iex,
    35		float *c1, float *c2, float *ryn, float *rzn)
    36	{
    37		for (int i = imin; i <  imax; i++) {
    38		for (int j = jmin; j <= jmax; j++) {
    39		for (int k = kmin; k <= kmax; k++) {
    40			int64_t n = (i * ni) + (j * nj) + (k * nk) + n0;
    41			updateEx_(
    42				n, nj, nk,
    43				ex, hy, hz,
    44				c1[iex[n]], c2[iex[n]], ryn[j], rzn[k]);
    45		}
    46		}
    47		}
    48	}
    49	
    50	void updateEx()
    51	{
    52		if (GPU) {
    53			dim3 grid(
    54				CEIL(kMax - kMin + 1, updateBlock.x),
    55				CEIL(jMax - jMin + 1, updateBlock.y),
    56				CEIL(iMax - iMin + 0, updateBlock.z));
    57			updateEx_gpu<<<grid, updateBlock>>>(
    58				iMin, jMin, kMin, iMax, jMax, kMax, Ni, Nj, Nk, N0,
    59				Ex, Hy, Hz, iEx,
    60				C1, C2, RYn, RZn);
    61			if (UM) cudaDeviceSynchronize();
    62		}
    63		else {
    64			updateEx_cpu(
    65				iMin, jMin, kMin, iMax, jMax, kMax, Ni, Nj, Nk, N0,
    66				Ex, Hy, Hz, iEx,
    67				C1, C2, RYn, RZn);
    68		}
    69	}

4.4.3 Mur吸収境界条件

Mur吸収境界条件はリスト4-4-2のようになります。
35,51,67行目のように1次元配列を1次元のblockとgridに分割しています。


リスト4-4-2 Mur吸収境界条件(murH.cu)
     1	__host__ __device__
     2	static void murh(float *h, mur_t *q, int64_t ni, int64_t nj, int64_t nk, int64_t n0)
     3	{
     4		int64_t m0 = (ni * q->i)  + (nj * q->j)  + (nk * q->k)  + n0;
     5		int64_t m1 = (ni * q->i1) + (nj * q->j1) + (nk * q->k1) + n0;
     6	
     7		h[m0] = q->f + q->g * (h[m1] - h[m0]);
     8		q->f = h[m1];
     9	}
    10	
    11	__global__
    12	static void murh_gpu(
    13		int64_t num, float *h, mur_t *fmur,
    14		int64_t ni, int64_t nj, int64_t nk, int64_t n0)
    15	{
    16		int64_t n = threadIdx.x + (blockIdx.x * blockDim.x);
    17		if (n < num) {
    18			murh(h, &fmur[n], ni, nj, nk, n0);
    19		}
    20	}
    21	
    22	static void murh_cpu(
    23		int64_t num, float *h, mur_t *fmur,
    24		int64_t ni, int64_t nj, int64_t nk, int64_t n0)
    25	{
    26		for (int64_t n = 0; n < num; n++) {
    27			murh(h, &fmur[n], ni, nj, nk, n0);
    28		}
    29	}
    30	
    31	void murHx()
    32	{
    33		assert(numMurHx > 0);
    34		if (GPU) {
    35			murh_gpu<<<(int)CEIL(numMurHx, murBlock), murBlock>>>(
    36				numMurHx, Hx, d_fMurHx,
    37				Ni, Nj, Nk, N0);
    38			if (UM) cudaDeviceSynchronize();
    39		}
    40		else {
    41			murh_cpu(
    42				numMurHx, Hx, fMurHx,
    43				Ni, Nj, Nk, N0);
    44		}
    45	}
    46	
    47	void murHy()
    48	{
    49		assert(numMurHy > 0);
    50		if (GPU) {
    51			murh_gpu<<<(int)CEIL(numMurHy, murBlock), murBlock>>>(
    52				numMurHy, Hy, d_fMurHy,
    53				Ni, Nj, Nk, N0);
    54			if (UM) cudaDeviceSynchronize();
    55		}
    56		else {
    57			murh_cpu(
    58				numMurHy, Hy, fMurHy,
    59				Ni, Nj, Nk, N0);
    60		}
    61	}
    62	
    63	void murHz()
    64	{
    65		assert(numMurHz > 0);
    66		if (GPU) {
    67			murh_gpu<<<(int)CEIL(numMurHz, murBlock), murBlock>>>(
    68				numMurHz, Hz, d_fMurHz,
    69				Ni, Nj, Nk, N0);
    70			if (UM) cudaDeviceSynchronize();
    71		}
    72		else {
    73			murh_cpu(
    74				numMurHz, Hz, fMurHz,
    75				Ni, Nj, Nk, N0);
    76		}
    77	}

4.4.4 PML吸収境界条件

PML吸収境界条件の内、例えばEx成分についてはリスト4-4-3のようになります。
60行目のように1次元配列を1次元のblockとgridに分割しています。


リスト4-4-3 PML吸収境界条件(Ex成分, pmlEx.cu)
     1	__host__ __device__
     2	static void pmlex(
     3		int64_t n, int64_t nj, int64_t nk,
     4		float *ex, float *hy, float *hz, float *exy, float *exz, float4 f)
     5	{
     6		*exy = (*exy + (f.x * (hz[n] - hz[n - nj]))) * f.z;
     7		*exz = (*exz - (f.y * (hy[n] - hy[n - nk]))) * f.w;
     8		ex[n] = *exy + *exz;
     9	}
    10	
    11	__global__
    12	static void pmlex_gpu(
    13		int ny, int nz,
    14		int64_t ni, int64_t nj, int64_t nk, int64_t n0,
    15		float *ex, float *hy, float *hz, float *exy, float *exz,
    16		int l, int64_t numpmlex,
    17		pml_t *fpmlex, float *rpmle, float *ryn, float *rzn, float *gpmlyn, float *gpmlzn)
    18	{
    19		int64_t n = threadIdx.x + (blockIdx.x * blockDim.x);
    20		if (n < numpmlex) {
    21			int i = fpmlex[n].i;
    22			int j = fpmlex[n].j;
    23			int k = fpmlex[n].k;
    24			int m = fpmlex[n].m;
    25			float ry = (m != PEC) ? ryn[MIN(MAX(j, 0), ny    )] * rpmle[m] : 0;
    26			float rz = (m != PEC) ? rzn[MIN(MAX(k, 0), nz    )] * rpmle[m] : 0;
    27			int64_t nc = (ni * i) + (nj * j) + (nk * k) + n0;
    28			pmlex(
    29				nc, nj, nk,
    30				ex, hy, hz, &exy[n], &exz[n],
    31				make_float4(ry, rz, gpmlyn[j + l], gpmlzn[k + l]));
    32		}
    33	}
    34	
    35	static void pmlex_cpu(
    36		int ny, int nz,
    37		int64_t ni, int64_t nj, int64_t nk, int64_t n0,
    38		float *ex, float *hy, float *hz, float *exy, float *exz,
    39		int l, int64_t numpmlex,
    40		pml_t *fpmlex, float *rpmle, float *ryn, float *rzn, float *gpmlyn, float *gpmlzn)
    41	{
    42		for (int64_t n = 0; n < numpmlex; n++) {
    43			int i = fpmlex[n].i;
    44			int j = fpmlex[n].j;
    45			int k = fpmlex[n].k;
    46			int m = fpmlex[n].m;
    47			float ry = (m != PEC) ? ryn[MIN(MAX(j, 0), ny    )] * rpmle[m] : 0;
    48			float rz = (m != PEC) ? rzn[MIN(MAX(k, 0), nz    )] * rpmle[m] : 0;
    49			int64_t nc = (ni * i) + (nj * j) + (nk * k) + n0;
    50			pmlex(
    51				nc, nj, nk,
    52				ex, hy, hz, &exy[n], &exz[n],
    53				make_float4(ry, rz, gpmlyn[j + l], gpmlzn[k + l]));
    54		}
    55	}
    56	
    57	void pmlEx()
    58	{
    59		if (GPU) {
    60			pmlex_gpu<<<(int)CEIL(numPmlEx, pmlBlock), pmlBlock>>>(
    61				Ny, Nz,
    62				Ni, Nj, Nk, N0,
    63				Ex, Hy, Hz, Exy, Exz,
    64				cPML.l, numPmlEx,
    65				d_fPmlEx, d_rPmlE, d_RYn, d_RZn, d_gPmlYn, d_gPmlZn);
    66			if (UM) cudaDeviceSynchronize();
    67		}
    68		else {
    69			pmlex_cpu(
    70				Ny, Nz,
    71				Ni, Nj, Nk, N0,
    72				Ex, Hy, Hz, Exy, Exz,
    73				cPML.l, numPmlEx,
    74				fPmlEx, rPmlE, RYn, RZn, gPmlYn, gPmlZn);
    75		}
    76	}

4.4.5 平均電磁界

リスト4-4-4に平均電磁界を計算する関数を示します。

ソースコードの説明
2行目:関数"average_e"は指定したセルの中心の電界3成分の絶対値の和を返します。
23行目:関数"average_h"は指定したセルの中心の磁界3成分の絶対値の和を返します。
46-48行目:blockを(Z,Y)方向の2次元、gridを(Z,Y,X)方向の3次元としています。
50行目:blockの大きさ(=スレッド数)を求めます。
51行目:スレッド番号を求めます。
52行目:ブロック番号を求めます。
60,68行目:スレッドの和(全体から見たら部分和)をreductionで求めます。
94-96行目:shared memoryは全スレッドで共有されるメモリーです。 reductionで必要になります。 shared memoryの大きさをexecution configulationの第3引数で渡します。 このときGPU側では44行目のように"extern __shared__"を付けます。
103-109行目:CPUで全体の和(=部分和の和)を計算します。


リスト4-4-4 平均電磁界(average.cu)
     1	__host__ __device__
     2	static float average_e(float *ex, float *ey, float *ez, int i, int j, int k, param_t *p)
     3	{
     4		return
     5			fabsf(
     6				+ ex[LA(p, i,     j,     k    )]
     7				+ ex[LA(p, i,     j + 1, k    )]
     8				+ ex[LA(p, i,     j,     k + 1)]
     9				+ ex[LA(p, i,     j + 1, k + 1)]) +
    10			fabsf(
    11				+ ey[LA(p, i,     j,     k    )]
    12				+ ey[LA(p, i,     j,     k + 1)]
    13				+ ey[LA(p, i + 1, j,     k    )]
    14				+ ey[LA(p, i + 1, j,     k + 1)]) +
    15			fabsf(
    16				+ ez[LA(p, i,     j,     k    )]
    17				+ ez[LA(p, i + 1, j,     k    )]
    18				+ ez[LA(p, i,     j + 1, k    )]
    19				+ ez[LA(p, i + 1, j + 1, k    )]);
    20	}
    21	
    22	__host__ __device__
    23	static float average_h(float *hx, float *hy, float *hz, int i, int j, int k, param_t *p)
    24	{
    25		return
    26			fabsf(
    27				+ hx[LA(p, i,     j,     k    )]
    28				+ hx[LA(p, i + 1, j,     k    )]) +
    29			fabsf(
    30				+ hy[LA(p, i,     j,     k    )]
    31				+ hy[LA(p, i,     j + 1, k    )]) +
    32			fabsf(
    33				+ hz[LA(p, i,     j,     k    )]
    34				+ hz[LA(p, i,     j,     k + 1)]);
    35	}
    36	
    37	// GPU
    38	__global__
    39	static void average_gpu(
    40		int nx, int ny, int nz, int imin, int imax,
    41		float *ex, float *ey, float *ez, float *hx, float *hy, float *hz,
    42		float *sume, float *sumh)
    43	{
    44		extern __shared__ float sm[];
    45	
    46		int i = imin + blockIdx.z;
    47		int j = threadIdx.y + (blockIdx.y * blockDim.y);
    48		int k = threadIdx.x + (blockIdx.x * blockDim.x);
    49	
    50		int nthread = blockDim.x * blockDim.y;
    51		int tid = threadIdx.x + (threadIdx.y * blockDim.x);
    52		int bid = blockIdx.x + (blockIdx.y * gridDim.x) + (blockIdx.z * gridDim.x * gridDim.y);
    53	
    54		if ((i < imax) && (j < ny) && (k < nz)) {
    55			sm[tid] = average_e(ex, ey, ez, i, j, k, &d_Param);
    56		}
    57		else {
    58			sm[tid] = 0;
    59		}
    60		reduction_sum(tid, nthread, sm, &sume[bid]);
    61	
    62		if ((i < imax) && (j < ny) && (k < nz)) {
    63			sm[tid] = average_h(hx, hy, hz, i, j, k, &d_Param);
    64		}
    65		else {
    66			sm[tid] = 0;
    67		}
    68		reduction_sum(tid, nthread, sm, &sumh[bid]);
    69	}
    70	
    71	// CPU
    72	static void average_cpu(float *sum)
    73	{
    74		sum[0] = 0;
    75		sum[1] = 0;
    76		for (int i = iMin; i < iMax; i++) {
    77		for (int j = 0; j < Ny; j++) {
    78		for (int k = 0; k < Nz; k++) {
    79			sum[0] += average_e(Ex, Ey, Ez, i, j, k, &h_Param);
    80			sum[1] += average_h(Hx, Hy, Hz, i, j, k, &h_Param);
    81		}
    82		}
    83		}
    84	}
    85	
    86	void average(double fsum[])
    87	{
    88		float sum[2];
    89	
    90		// sum
    91		if (GPU) {
    92			cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
    93	
    94			int sm_size = sumBlock.x * sumBlock.y * sizeof(float);
    95			average_gpu<<<sumGrid, sumBlock, sm_size>>>(
    96				Nx, Ny, Nz, iMin, iMax, Ex, Ey, Ez, Hx, Hy, Hz, m_sumE, m_sumH);
    97
    98			// device to host
    99			size_t size = sumGrid.x * sumGrid.y * sumGrid.z * sizeof(float);
   100			cudaMemcpy(h_sumE, d_sumE, size, cudaMemcpyDeviceToHost);
   101			cudaMemcpy(h_sumH, d_sumH, size, cudaMemcpyDeviceToHost);
   102	
   103			// partial sum
   104			sum[0] = 0;
   105			sum[1] = 0;
   106			for (int i = 0; i < size / sizeof(float); i++) {
   107				sum[0] += m_sumE[i];
   108				sum[1] += m_sumH[i];
   109			}
   110		}
   111		else {
   112			average_cpu(sum);
   113		}
   114	
   115		// average
   116		fsum[0] = sum[0] / (4.0 * Nx * Ny * Nz);
   117		fsum[1] = sum[1] / (2.0 * Nx * Ny * Nz);
   118	}

リスト4-4-5にreduction関数を示します。
reductionとはN個の演算を複数のスレッドの共同作業でlog2N回で求めるGPUの操作です。
引数"s"がshared memoryでありスレッド間で共有されます。
図4-4-1にreductionの概念図を示します。
詳細についてはCUDAのサンプルプログラム"reduction"を参考にして下さい。


リスト4-4-5 reduction関数
     1	__device__
     2	void reduction_sum(int tid, int n, float *s, float *sum)
     3	{
     4		__syncthreads();
     5	
     6		for (int stride = (n + 1) >> 1; n > 1; stride = (stride + 1) >> 1, n = (n + 1) >> 1) {
     7			if (tid + stride < n) {
     8				s[tid] += s[tid + stride];
     9			}
    10			__syncthreads();
    11		}
    12	
    13		if (tid == 0) {
    14			*sum = s[0];
    15		}
    16	}


図4-4-1 reductionの概念図

4.4.6 CUDAの計算時間

表4-4-1にCPU(40コア)とGPUの計算時間を示します。
GPUはCPU40コアの4〜5倍速くなり、問題のサイズが大きくなるほど速度比も上がります。

表4-4-1 CUDAの計算時間
(FOCUSスパコンFシステム、()内はCPU40コアとの速度比)
CPU/GPUベンチマーク200ベンチマーク300ベンチマーク400ベンチマーク500
CPU+MPI(40コア)19.3秒 (1.0)71.1秒 (1.0)181.4秒 (1.0)402.0秒 (1.0)
GPU5.4秒 (3.6)17.2秒 (4.1)39.9秒 (4.5)79.4秒 (5.1)

4.4.7 ブロックの大きさ

リスト4-4-1の57-59行目のブロックの大きさによって計算時間が変わり、 ブロックの大きさには最適値が存在します。
図4-4-2にblockとgridの概念図を示します。
表4-4-2にblockDim.xとblockDim.yを変えたときの計算時間を示します。
ブロックサイズ(= blockDim.x * blockDim.y)が32以上であれば計算時間に大きな違いはありませんが、 11秒以下の組み合わせに色をつけた通り、ブロックサイズが64〜256のとき計算時間は最小になります。 この最適値は問題のサイズにはあまり関係しません。 なお、ブロックサイズの上限はCompute Capability 2.0以降では1024です。
以上の結果から本プログラムでは blockDim.x = 32, blockDim.y = 4 とします。 前者が最も小さい組み合わせを選んだ理由は、Z方向のセル数Nzが常に十分大きいとは限らないので、 Nzが小さいときはblockDim.xが小さい方が効率よく計算できるためです。


図4-4-2 FDTD法におけるblockとgridの概念図

表4-4-2 CUDAの計算時間とblockの関係
(Windows、GTX1070、ベンチマーク200、blockDim.z=1)
blockDim.yblockDim.x
163264128256
115.0秒11.2秒10.9秒10.6秒10.6秒
211.3秒10.7秒10.7秒10.8秒10.8秒
411.1秒10.7秒10.8秒11.1秒12.0秒
811.3秒10.8秒11.1秒11.9秒N/A
1611.4秒11.1秒11.8秒N/AN/A

4.4.8 プロファイル

コマンドラインで以下を実行するとプロファイルが標準エラーに出力されます。
$ nvprof ofd_cuda データファイル
表4-4-3にプロファイルの一例を示します。 電磁界6成分の更新が計算時間の96%を占めることがわかります。

表4-4-3 GPUプログラムのプロファイル
            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:   16.26%  6.51108s      2001  3.2539ms  3.2263ms  3.3336ms  updateEz_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                   16.19%  6.48402s      2001  3.2404ms  3.2087ms  3.3276ms  updateEx_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                   16.18%  6.47927s      2001  3.2380ms  3.2111ms  3.3234ms  updateHy_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                   15.93%  6.37880s      2001  3.1878ms  3.1378ms  3.2765ms  updateEy_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                   15.90%  6.36630s      2001  3.1816ms  3.1409ms  3.2716ms  updateHz_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                   15.80%  6.32781s      2001  3.1623ms  3.1214ms  3.2530ms  updateHx_gpu(int, int, int, int, int, int, long, long, long, long, float*, float*, float*, unsigned char*, float*, float*, float*, float*)
                    3.30%  1.32047s      6003  219.97us  141.73us  271.45us  murh_gpu(long, float*, mur_t*, long, long, long, long)
                    0.26%  102.75ms        11  9.3408ms  9.2508ms  10.230ms  average_gpu(int, int, int, int, int, float*, float*, float*, float*, float*, float*, float*, float*)
                    0.14%  57.293ms      2032  28.195us     832ns  9.3392ms  [CUDA memcpy HtoD]
                    0.04%  14.926ms      2001  7.4590us  5.9840us  8.5120us  efeed_gpu(int, feed_t*, float*, float*, float, float, int, int, float, float, float*, float*, float*, float*, float*, float*)
                    0.00%  1.7937ms        24  74.736us  1.5680us  83.583us  [CUDA memcpy DtoH]
                    0.00%  41.343us        40  1.0330us     800ns  1.5680us  [CUDA memset]
      API calls:   98.27%  39.8028s      2012  19.783ms  9.3660us  22.910ms  cudaMemcpyToSymbol
                    0.84%  338.46ms        30  11.282ms  5.1380us  330.43ms  cudaMalloc
                    0.41%  165.22ms        44  3.7551ms  7.0590us  10.907ms  cudaMemcpy
                    0.32%  128.87ms     20021  6.4360us  4.9970us  444.91us  cudaLaunchKernel
                    0.16%  64.031ms        28  2.2868ms  4.1970us  6.5694ms  cudaFree
                    0.00%  749.90us        40  18.747us  2.7750us  75.324us  cudaMemset