図4-5-1 reductionの概念図
CUDAによるFDTD法の高速化では下記の点に注意してプログラムを作成します。
リスト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では最後に同期を取っています。
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-5-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-5-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-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で全体の和(=部分和の和)を計算します。
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"を参考にして下さい。
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にWindowsでの計算時間の測定結果を示します。
参考までに表4-2-1のCPU1コア(SIMDなし)の結果も載せます。
GPUを用いるとCPU1コア(SIMDなし)の11〜29倍速くなることがわかります。
一般に問題のサイズが大きくなるほど速度比も上がります。
| ベンチマークNo. | CPU(1コア,SIMDなし) | GPU |
|---|---|---|
| 1 | 29.1秒(1.0) | 2.6秒(11.2) |
| 2 | 233.5秒(1.0) | 10.8秒(21.6) |
| 3 | 2476.3秒(1.0) | 84.4秒(29.3) |
リスト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が小さい方が効率よく計算できるためです。
| 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 |
図4-5-3にNVIDIA Visual Profilerの出力図の一例を示します。
電磁界6成分の更新がほぼすべての計算時間を占め、
その間のアイドルがなく効率よく計算されていることがわかります。
