4.5 CUDAによる高速化

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

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

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

4.5.2 電磁界の更新

リスト4-5-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-5-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.5.3 Mur吸収境界条件

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


リスト4-5-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.5.4 PML吸収境界条件

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


リスト4-5-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.5.5 平均電磁界

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


リスト4-5-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-5-1 reductionの概念図

4.5.6 CUDAの計算時間

表4-5-1にWindowsでの計算時間の測定結果を示します。
参考までに表4-2-1のCPU1コア(SIMDなし)の結果も載せます。
GPUを用いるとCPU1コア(SIMDなし)の11〜29倍速くなることがわかります。 一般に問題のサイズが大きくなるほど速度比も上がります。

表4-5-1 CUDAの計算時間(Windows, GTX1070(1920コア、GDDR5 8GB)、()内はCPU1コアとの速度比)
ベンチマークNo.CPU(1コア,SIMDなし)GPU
129.1秒(1.0)2.6秒(11.2)
2233.5秒(1.0)10.8秒(21.6)
32476.3秒(1.0)84.4秒(29.3)

4.5.7 ブロックの大きさ

リスト4-5-1の57-59行目のブロックの大きさによって計算時間が変わり、 ブロックの大きさには最適値が存在します。
図4-5-2にblockとgridの概念図を示します。
表4-5-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-5-2 FDTD法におけるblockとgridの概念図

表4-5-2 CUDAの計算時間とblockの関係(Windows、GTX1070、ベンチマークNo.2、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.5.8 プロファイル

図4-5-3にNVIDIA Visual Profilerの出力図の一例を示します。
電磁界6成分の更新がほぼすべての計算時間を占め、 その間のアイドルがなく効率よく計算されていることがわかります。


図4-5-3 NVIDIA Visual Profilerの出力図(ベンチマークNo.2、GTX1070)