3.5 CUDAによる並列化

3.5.1 行列・ベクトル積

行列とベクトルの積を計算するには、(i,j,k)の三重ループをgrid,block分解して行います。 その際、blockの内側ループがアクセスするアドレスが連続して格納されていることが重要です。

3.5.2 ベクトル演算

CUDAによるベクトルのコピーはリスト3-5-1の通りです。

リスト3-5-1 CUDAによるコピー(Zcopy)

__global__
static void Zcopy_kernel(int64_t n, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	int64_t tid = threadIdx.x + (blockIdx.x * blockDim.x);

	if (tid < n) {
		y[tid].x = x[tid].x;
		y[tid].y = x[tid].y;
	}
}

void Zcopy(int64_t n, const cuDoubleComplex *x, cuDoubleComplex *y)
{
	unsigned int grid = (unsigned int)CEIL(n, BLAS_BLOCK_SIZE);
	Zcopy_kernel<<<grid, BLAS_BLOCK_SIZE>>>(n, x, y);
}

CUDAによる内積の計算はリスト3-5-2の通りです。
ここでreduction_sumはshared memoryに格納された配列の和をredunction演算によって求める関数です。
また、配列Wblasは部分和を格納する作業配列です。

リスト3-5-2 CUDAによる内積(Zdotu)

__global__
static void Zdotu_kernel(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y, cuDoubleComplex *work)
{
	extern __shared__ double s[];

	int64_t tid = (blockIdx.x * blockDim.x) + threadIdx.x;

	double sum_r = 0;
	double sum_i = 0;
	while (tid < n) {
		sum_r += (x[tid].x * y[tid].x) - (x[tid].y * y[tid].y);
		sum_i += (x[tid].x * y[tid].y) + (x[tid].y * y[tid].x);
		tid += blockDim.x * gridDim.x;
	}
	s[0          + threadIdx.x] = sum_r;
	s[blockDim.x + threadIdx.x] = sum_i;

	reduction_sum(threadIdx.x, blockDim.x, &s[0         ], &work[blockIdx.x].x);
	reduction_sum(threadIdx.x, blockDim.x, &s[blockDim.x], &work[blockIdx.x].y);
}

cuDoubleComplex Zdotu(int64_t n, const cuDoubleComplex *x, const cuDoubleComplex *y)
{
	cuDoubleComplex ret = make_cuDoubleComplex(0, 0);

	size_t sm_size = 2 * BLAS_BLOCK_SIZE * sizeof(double);
	Zdotu_kernel<<<BLAS_GRID_SIZE, BLAS_BLOCK_SIZE, sm_size>>>(n, x, y, d_Wblas);
	cudaMemcpy(h_Wblas, d_Wblas, BLAS_GRID_SIZE * sizeof(cuDoubleComplex), cudaMemcpyDeviceToHost);
	for (int i = 0; i < BLAS_GRID_SIZE; i++) {
		ret = cuCadd(ret, h_Wblas[i]);
	}

	return ret;
}

3.5.3 GPUプロファイル

Windows環境では図3-5-1のようにNVIDIA Visual Profilerを用いるとプロファイルが得られます。 右側は反復回数1回のtimelineです。

図3-5-1 NVIDIA Visual Profilerの出力図(nomatrixモード)

WindowsとLinux環境ではコマンドラインツールnvprofを用いて

$ nvprof oth_cuda [-maxtix|-nomatrix] 入力データ

と行います。図3-5-2に出力結果を示します。
これからmatrixモードでは行列・ベクトル積に64%、ベクトル演算に35%、 nomatrixモードでは行列・ベクトル積に68%、ベクトル演算に32%の時間が消費されていることがわかります。

            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:   21.53%  1.33839s      1003  1.3344ms  1.3141ms  1.3566ms  prodmvEz_gpu
                   21.38%  1.32862s      1003  1.3246ms  1.3023ms  1.3501ms  prodmvEx_gpu
                   21.37%  1.32843s      1003  1.3245ms  1.3058ms  1.3447ms  prodmvEy_gpu
                   11.11%  690.52ms      2512  274.89us  172.32us  481.63us  Zdotc_kernel
                   10.71%  665.96ms      1001  665.29us  662.26us  673.21us  Zaxbycz_kernel
                    8.20%  509.79ms      1003  508.27us  501.95us  520.38us  Zaxpy_kernel
                    5.46%  339.58ms      1005  337.89us  331.84us  348.54us  Zcopy_kernel

(a) matrixモード

            Type  Time(%)      Time     Calls       Avg       Min       Max  Name
 GPU activities:   22.93%  1.58945s      1004  1.5831ms  1.1302ms  1.8267ms  matrixEy_gpu
                   22.65%  1.56996s      1004  1.5637ms  1.1388ms  1.8295ms  matrixEx_gpu
                   22.44%  1.55543s      1004  1.5492ms  1.0794ms  1.7497ms  matrixEz_gpu
                   10.01%  694.19ms      2512  276.35us  172.38us  360.86us  Zdotc_kernel
                    9.61%  666.11ms      1001  665.45us  662.26us  672.98us  Zaxbycz_kernel
                    7.35%  509.65ms      1003  508.12us  501.63us  519.35us  Zaxpy_kernel
                    4.90%  339.61ms      1005  337.92us  336.16us  348.86us  Zcopy_kernel

(b) nomatrixモード
図3-5-2 nvprofの出力