図4-5-1 reductionの概念図
CUDAによるFDTD法の高速化では下記の点に注意してプログラムを作成します。
ここではCUDA6以降のUnified Memoryについて説明します。
リスト4-5-1にEx成分を更新する関数を示します。
[3]の手法を用いてCPUコードとGPUコードを共通化しています。
ソースコードの説明
1行目:計算の実体である関数はCPUとGPUで共通です。
従って関数修飾子"__host__ __device__"を付けます。
12行目:GPUコードには関数修飾子"__global__"を付けます。
19-21行目:gridとblockからi,j,kを計算します。
23-25行目:i,j,kのオーバーフローをチェックします(セル数がブロックサイズの倍数とは限らないため)。
56-58行目:blockは内側から順にk,j,iとします。
アドレスはkについて連続ですのでこのようにする必要があります(coalescing of memory)。
1 __host__ __device__
2 static void updateEx_(
3 int64_t n, int64_t nj, int64_t 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 nx, int ny, int nz,
15 int64_t ni, int64_t nj, int64_t nk, int64_t n0,
16 float *ex, float *hy, float *hz, unsigned char *iex,
17 float *c1, float *c2, float *ryn, float *rzn)
18 {
19 int i = threadIdx.z + (blockIdx.z * blockDim.z);
20 int j = threadIdx.y + (blockIdx.y * blockDim.y);
21 int k = threadIdx.x + (blockIdx.x * blockDim.x);
22
23 if ((i < nx) &&
24 (j <= ny) &&
25 (k <= nz)) {
26 int64_t n = (i * ni) + (j * nj) + (k * nk) + n0;
27 updateEx_(
28 n, nj, nk,
29 ex, hy, hz,
30 c1[iex[n]], c2[iex[n]], ryn[j], rzn[k]);
31 }
32 }
33
34 static void updateEx_cpu(
35 int nx, int ny, int nz,
36 int64_t ni, int64_t nj, int64_t nk, int64_t n0,
37 float *ex, float *hy, float *hz, unsigned char *iex,
38 float *c1, float *c2, float *ryn, float *rzn)
39 {
40 for (int i = 0; i < nx; i++) {
41 for (int j = 0; j <= ny; j++) {
42 for (int k = 0; k <= nz; k++) {
43 int64_t n = (i * ni) + (j * nj) + (k * nk) + n0;
44 updateEx_(
45 n, nj, nk,
46 ex, hy, hz,
47 c1[iex[n]], c2[iex[n]], ryn[j], rzn[k]);
48 }
49 }
50 }
51 }
52
53 void updateEx()
54 {
55 if (GPU) {
56 dim3 grid(((Nz + 1) + (updateBlock.x - 1)) / updateBlock.x,
57 ((Ny + 1) + (updateBlock.y - 1)) / updateBlock.y,
58 ((Nx + 0) + (updateBlock.z - 1)) / updateBlock.z);
59 updateEx_gpu<<<grid, updateBlock>>>(
60 Nx, Ny, Nz,
61 Ni, Nj, Nk, N0,
62 Ex, Hy, Hz, iEx,
63 C1, C2, RYn, RZn);
64 cudaDeviceSynchronize();
65 }
66 else {
67 updateEx_cpu(
68 Nx, Ny, Nz,
69 Ni, Nj, Nk, N0,
70 Ex, Hy, Hz, iEx,
71 C1, C2, RYn, RZn);
72 }
73 }
Mur吸収境界条件はリスト4-5-2のようになります。
32,46,60行目のように1次元配列を1次元のblockとgridに分割しています。
1 __host__ __device__
2 static void murh(float *h, mur2_t *q, param_t *p)
3 {
4 int64_t n0 = LA(p, q->i, q->j, q->k );
5 int64_t n1 = LA(p, q->i1, q->j1, q->k1);
6
7 h[n0] = q->f + q->g * (h[n1] - h[n0]);
8 q->f = h[n1];
9 }
10
11 __global__
12 static void murh_gpu(int64_t num, float *h, mur2_t *fmur)
13 {
14 int64_t n = threadIdx.x + (blockIdx.x * blockDim.x);
15 if (n < num) {
16 murh(h, &fmur[n], &d_Param);
17 }
18 }
19
20 static void murh_cpu(int64_t num, float *h, mur2_t *fmur)
21 {
22 for (int64_t n = 0; n < num; n++) {
23 murh(h, &fmur[n], &h_Param);
24 }
25 }
26
27 void murHx_gpu()
28 {
29 assert(numMurHx > 0);
30 if (GPU) {
31 cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
32 int grid = (numMurHx + (murBlock - 1)) / murBlock;
33 murh_gpu<<<grid, murBlock>>>(numMurHx, Hx, fMurHx);
34 cudaDeviceSynchronize();
35 }
36 else {
37 murh_cpu(numMurHx, Hx, fMurHx);
38 }
39 }
40
41 void murHy_gpu()
42 {
43 assert(numMurHy > 0);
44 if (GPU) {
45 cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
46 int grid = (numMurHy + (murBlock - 1)) / murBlock;
47 murh_gpu<<<grid, murBlock>>>(numMurHy, Hy, fMurHy);
48 cudaDeviceSynchronize();
49 }
50 else {
51 murh_cpu(numMurHy, Hy, fMurHy);
52 }
53 }
54
55 void murHz_gpu()
56 {
57 assert(numMurHz > 0);
58 if (GPU) {
59 cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
60 int grid = (numMurHz + (murBlock - 1)) / murBlock;
61 murh_gpu<<<grid, murBlock>>>(numMurHz, Hz, fMurHz);
62 cudaDeviceSynchronize();
63 }
64 else {
65 murh_cpu(numMurHz, Hz, fMurHz);
66 }
67 }
PML吸収境界条件の内、例えばEx成分についてはリスト4-5-3のようになります。
62行目のように1次元配列を1次元のblockとgridに分割しています。
1 __host__ __device__
2 static void pmlex(
3 int i, int j, int k,
4 float *ex, float *hy, float *hz, float *exy, float *exz,
5 float rm1, float rm2, float gyn, float gzn, param_t *p)
6 {
7 float dhz = hz[LA(p, i, j, k)] - hz[LA(p, i, j - 1, k)];
8 *exy = (*exy + (rm1 * dhz)) * gyn;
9
10 float dhy = hy[LA(p, i, j, k)] - hy[LA(p, i, j, k - 1)];
11 *exz = (*exz - (rm2 * dhy)) * gzn;
12
13 ex[LA(p, i, j, k)] = *exy + *exz;
14 }
15
16 __global__
17 static void pmlex_gpu(
18 int ny, int nz,
19 float *ex, float *hy, float *hz, float *exy, float *exz,
20 int64_t numpmlex, pml2_t *fpmlex, float *rpmle, float *ryn, float *rzn,
21 float *gpmlyn, float *gpmlzn, int l)
22 {
23 int n64_t = threadIdx.x + (blockIdx.x * blockDim.x);
24 if (n < numpmlex) {
25 int i = fpmlex[n].i;
26 int j = fpmlex[n].j;
27 int k = fpmlex[n].k;
28 int m = fpmlex[n].m;
29 float rm1 = (m != PEC) ? ryn[MIN(MAX(j, 0), ny )] * rpmle[m] : 0;
30 float rm2 = (m != PEC) ? rzn[MIN(MAX(k, 0), nz )] * rpmle[m] : 0;
31 pmlex(
32 i, j, k,
33 ex, hy, hz, &exy[n], &exz[n],
34 rm1, rm2, gpmlyn[j + l], gpmlzn[k + l], &d_Param);
35 }
36 }
37
38 static void pmlex_cpu(
39 int ny, int nz,
40 float *ex, float *hy, float *hz, float *exy, float *exz,
41 int64_t numpmlex, pml2_t *fpmlex, float *rpmle, float *ryn, float *rzn,
42 float *gpmlyn, float *gpmlzn, int l)
43 {
44 for (int64_t n = 0; n < numpmlex; n++) {
45 int i = fpmlex[n].i;
46 int j = fpmlex[n].j;
47 int k = fpmlex[n].k;
48 int m = fpmlex[n].m;
49 float rm1 = (m != PEC) ? ryn[MIN(MAX(j, 0), ny )] * rpmle[m] : 0;
50 float rm2 = (m != PEC) ? rzn[MIN(MAX(k, 0), nz )] * rpmle[m] : 0;
51 pmlex(
52 i, j, k,
53 ex, hy, hz, &exy[n], &exz[n],
54 rm1, rm2, gpmlyn[j + l], gpmlzn[k + l], &h_Param);
55 }
56 }
57
58 void pmlEx_gpu()
59 {
60 if (GPU) {
61 cudaMemcpyToSymbol(d_Param, &h_Param, sizeof(param_t));
62 int grid = (numPmlEx + (pmlBlock - 1)) / pmlBlock;
63 pmlex_gpu<<<grid, pmlBlock>>>(
64 Ny, Nz,
65 Ex, Hy, Hz, Exy, Exz,
66 numPmlEx, fPmlEx, rPmlE, RYn, RZn,
67 gPmlYn, gPmlZn, cPML.l);
68 cudaDeviceSynchronize();
69 }
70 else {
71 pmlex_cpu(
72 Ny, Nz,
73 Ex, Hy, Hz, Exy, Exz,
74 numPmlEx, fPmlEx, rPmlE, RYn, RZn,
75 gPmlYn, gPmlZn, cPML.l);
76 }
77 }
リスト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__"を付けます。
100-106行目: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,
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 = 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 < nx) && (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 < nx) && (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 = 0; i < Nx; 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, Ex, Ey, Ez, Hx, Hy, Hz, m_sumE, m_sumH);
97 cudaDeviceSynchronize();
98
99 // partial sum
100 sum[0] = 0;
101 sum[1] = 0;
102 int n = sumGrid.x * sumGrid.y * sumGrid.z;
103 for (int i = 0; i < n; i++) {
104 sum[0] += m_sumE[i];
105 sum[1] += m_sumH[i];
106 }
107 }
108 else {
109 average_cpu(sum);
110 }
111
112 // average
113 fsum[0] = sum[0] / (4.0 * Nx * Ny * Nz);
114 fsum[1] = sum[1] / (2.0 * Nx * Ny * Nz);
115 }
リスト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コアの8-11倍速くなることがわかります。
一般に問題のサイズが大きくなるほど速度比も上がります。
| ベンチマークNo. | CPU(1コア) | GPU |
|---|---|---|
| 1 | 27.6秒(1.0) | 3.3秒(8.4) |
| 2 | 229.6秒(1.0) | 21.4秒(10.7) |
リスト4-5-1の60-62行目のブロックの大きさによって計算時間が変わり、
ブロックの大きさには最適値が存在します。
図4-5-2にblockとgridの概念図を示します。
表4-5-2にblockDim.xとblockDim.yを変えたときの計算時間を示します。
色のついた組み合わせのとき計算時間が最小になります。
すなわちblockDim.xは32以上、
ブロックサイズ(= blockDim.x * blockDim.y)が128か256が最適値となります。
この最適値は問題のサイズにはあまり関係しません。
なお、ブロックサイズの上限はCompute Capability 2.0以降では1024です。
以上の結果から本プログラムでは blockDim.x = 32, blockDim.y = 8 とします。
前者が最も小さい組み合わせを選んだ理由は、Z方向のセル数Nzが常に十分大きいとは限らないので、
Nzが小さいときはblockDim.xが小さい方が効率よく計算できるためです。
| blockDim.y | blockDim.x | ||||
|---|---|---|---|---|---|
| 16 | 32 | 64 | 128 | 256 | |
| 1 | 76.8秒 | 43.9秒 | 28.5秒 | 22.3秒 | 22.4秒 |
| 2 | 44.2秒 | 26.0秒 | 22.4秒 | 22.9秒 | 25.3秒 |
| 4 | 28.6秒 | 22.2秒 | 22.9秒 | 25.6秒 | 29.5秒 |
| 8 | 29.0秒 | 22.5秒 | 25.8秒 | 31.8秒 | N/A |
| 16 | 29.0秒 | 24.1秒 | 31.6秒 | N/A | N/A |