図4-4-1 reductionの概念図
CUDAによるFDTD法の高速化では下記の点に注意してプログラムを作成します。
リスト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では最後に同期を取っています。
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 }
Mur吸収境界条件はリスト4-4-2のようになります。
35,51,67行目のように1次元配列を1次元のblockとgridに分割しています。
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 }
PML吸収境界条件の内、例えばEx成分についてはリスト4-4-3のようになります。
60行目のように1次元配列を1次元のblockとgridに分割しています。
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-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で全体の和(=部分和の和)を計算します。
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"を参考にして下さい。
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の57-59行目のブロックの大きさによって計算時間が変わり、
ブロックの大きさには最適値が存在します。
図4-4-2にblockとgridの概念図を示します。
表4-4-1に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が小さい方が効率よく計算できるためです。
| blockDim.y | blockDim.x | ||||
|---|---|---|---|---|---|
| 16 | 32 | 64 | 128 | 256 | |
| 1 | 15.0秒 | 11.2秒 | 10.9秒 | 10.6秒 | 10.6秒 |
| 2 | 11.3秒 | 10.7秒 | 10.7秒 | 10.8秒 | 10.8秒 |
| 4 | 11.1秒 | 10.7秒 | 10.8秒 | 11.1秒 | 12.0秒 |
| 8 | 11.3秒 | 10.8秒 | 11.1秒 | 11.9秒 | N/A |
| 16 | 11.4秒 | 11.1秒 | 11.8秒 | N/A | N/A |
コマンドラインで以下を実行するとプロファイルが標準エラーに出力されます。
$ nvprof ofd_cuda データファイル
表4-4-2にプロファイルの一例を示します。
電磁界6成分の更新が計算時間の96%を占めることがわかります。
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