3.4 GPUによる並列計算

3.4.1 GPUによる並列計算

リスト3-4-1にGPU版のred-black SOR法の反復計算部分を示します。
blockは(k, j)の2次元の(16, 16)とし、gridは((Nz+1)/16, (Ny+1)/16, Nx+1)の3次元とします。
これは(i, j, k)の3次元空間を1次元配列で並べるときに内側からk, j, iの順にとっているためです。 CUDAでは変数のアドレスがスレッド番号の順に並んでいることが重要です[8]。
残差の計算はblock内で部分和をshared memoryのreduction演算で求め、 部分和の配列をCPUに転送しその和をCPUで計算しています。 このようにすると変数全体を転送して残差を計算することに比べると転送量が1/256倍になり、 データ転送に要する時間を無視することができます。
GPU版をMPIと併用することにより複数のグラフィックボードを用いて並列計算することができます。


リスト3-4-1 CUDAを用いた反復計算部の並列化(sor.cuの一部, GPU版, novectorモード)
     1	#include "ost.h"
     2	#include "ost_cuda.h"
     3	#include "ost_prototype.h"
     4	#include "reduction_sum.cu"
     5	
     6	__host__ __device__
     7	static real_t update(
     8		int i, int j, int k, int64_t ni, int64_t nj, int64_t nk, int64_t n,
     9		real_t *rxp, real_t *ryp, real_t *rzp, real_t *rxm, real_t *rym, real_t *rzm,
    10		real_t *v, id_t *idvolt, id_t *idepsr, real_t *epsr, real_t omega)
    11	{
    12		const real_t e = epsr[idepsr[n]];
    13	
    14		const real_t axp = (e + epsr[idepsr[n + ni]]) * rxp[i];
    15		const real_t axm = (e + epsr[idepsr[n - ni]]) * rxm[i];
    16		const real_t ayp = (e + epsr[idepsr[n + nj]]) * ryp[j];
    17		const real_t aym = (e + epsr[idepsr[n - nj]]) * rym[j];
    18		const real_t azp = (e + epsr[idepsr[n + nk]]) * rzp[k];
    19		const real_t azm = (e + epsr[idepsr[n - nk]]) * rzm[k];
    20	
    21		const real_t res = (idvolt[n] ? 0 : 1) * omega * ((
    22			axp * v[n + ni] +
    23			axm * v[n - ni] +
    24			ayp * v[n + nj] +
    25			aym * v[n - nj] +
    26			azp * v[n + nk] +
    27			azm * v[n - nk]) / (axp + axm + ayp + aym + azp + azm) - v[n]);
    28	
    29		v[n] += res;
    30	
    31		return (res * res);
    32	}
    33	
    34	__global__
    35	static void update_gpu(int oe,
    36		int nx, int ny, int nz,
    37		int64_t ni, int64_t nj, int64_t nk,
    38		real_t *rxp, real_t *ryp, real_t *rzp, real_t *rxm, real_t *rym, real_t *rzm,
    39		real_t *v, id_t *idvolt, id_t *idepsr, real_t *epsr, real_t omega, real_t *d_res2)
    40	{
    41		extern __shared__ real_t sm[];
    42	
    43		const int i = blockIdx.z;
    44		const int j = threadIdx.y + (blockIdx.y * blockDim.y);
    45		const int k = threadIdx.x + (blockIdx.x * blockDim.x);
    46	
    47		const int nthread = blockDim.x * blockDim.y;
    48		const int tid = threadIdx.x + (threadIdx.y * blockDim.x);
    49		const int bid = blockIdx.x + (blockIdx.y * gridDim.x) + (blockIdx.z * gridDim.x * gridDim.y);
    50	
    51		real_t res2 = 0;
    52	
    53		if (i <= nx) {
    54		if (j <= ny) {
    55		if (k <= nz) {
    56			if ((i + j + k) % 2 == oe) {
    57				const int64_t n = ((i + 1) * ni) + ((j + 1) * nj) + ((k + 1) * nk);
    58				res2 = update(
    59					i, j, k, ni, nj, nk, n,
    60					rxp, ryp, rzp, rxm, rym, rzm,
    61					v, idvolt, idepsr, epsr, omega);
    62			}
    63		}
    64		}
    65		}
    66		sm[tid] = res2;
    67	
    68		reduction_sum(tid, nthread, sm, &d_res2[bid]);
    69	}


3.4.2 GPUの計算時間

表3-4-1にノード数を変えたときの計算時間を示します。
2ノードは1ノードの2倍近く速くなります。

表3-4-1 GPUの計算時間(FOCUSスパコンFシステム, benchmark600, novectorモード)
ノード数計算時間(速度比)
136.7秒 (1.0)
220.0秒 (1.84)