GPGPU ライブラリ “ArrayFire” を生 CUDA と比較する

後輩君が研究室のサーバーに ArrayFire を入れてくれた． GPGPU のライブラリは大抵「CPU と比べるととても速い」と謳っているが，実際は直に CUDA で組んだ場合と比べてどうなのかが一番重要なので，簡単な例で試してみることにした．

ArrayFire について

ArrayFire はいわゆる GPGPU 用の C/C++ ライブラリ．これまでしばらくの間商用ライセンスと一部フリーで提供されていたが， 2014年11月にオープンソースとして公開され，今は github から clone してインストールできる．バックグラウンドには CUDA のほか OpenCL も選べるので，MacOSX でも動かせるようだ．

ArrayFire の特徴は，GPU 上のデータ構造を単一の array というクラスで表しているところにある． array は GPU メモリ上のオブジェクトをラップしていて，メモリの確保や転送を隠ぺいしてくれる上に，dim を設定することで行列や多次元配列のように扱うこともできる．さらに，基本的な行列演算や数値計算を行う関数が用意されていて，これを使えば GPU の知識がなくても GPU 上での計算が実行できる．

で，GPGPU を使う側として気になるのは，メモリ転送やアルゴリズムを隠ぺいした分，CUDA でハンドチューニングした場合に比べてどれだけ速度が遅くなるかというところ．足し算引き算などの element-base な演算だとさすがに差がつかないと思うので，reduction の計算と，ストリームを使った非同期計算を試してみることにした．

Reduction (sum) のパフォーマンス

Reduction の計算として，単純なベクトル要素の総和の計算を行ってみた．

生 CUDA では，共有メモリを使った一般的な実装を行った．

https://github.com/hnagata/af-performance/blob/master/cuda_sum.cu

#include <vector>
#include <cuda.h>
#include "perform.h"
#include "cuda_common.h"

const int GRID_SIZE = 448;
const int BLOCK_SIZE = 32;

__global__ void sum_dev(double *a, int n, double *out) {
	__shared__ double work[BLOCK_SIZE];
	const int tidx = threadIdx.x;
	const int bidx = blockIdx.x;
	const int idx = threadIdx.x + blockIdx.x * blockDim.x;
	const int dim = blockDim.x * gridDim.x;

	work[tidx] = 0;

	for (int i = idx; i < n; i += dim) {
		work&#91;tidx&#93; += a&#91;i&#93;;
	}
	__syncthreads();

	int w = BLOCK_SIZE;
	while (w > 2) {
		w /= 2;
		if (tidx < w) work&#91;tidx&#93; = work&#91;tidx + w&#93;;
		__syncthreads();
	}

	if (tidx == 0) {
		out&#91;bidx&#93; = work&#91;0&#93; + work&#91;1&#93;;
	}
}

double sum(const std::vector<double> &a) {
	double *a_dev = NULL;
	double ans = 0;	
	CUDA_CALL(cudaMalloc((void**) &a_dev, sizeof(double) * a.size()));
	CUDA_CALL(cudaMemcpy(
			a_dev, a.data(), sizeof(double) * a.size(), 
			cudaMemcpyHostToDevice));
	sum_dev<<<GRID_SIZE, BLOCK_SIZE>>>(a_dev, a.size(), a_dev);
	CUDA_CHECK();
	sum_dev<<<1, BLOCK_SIZE>>>(a_dev, GRID_SIZE, a_dev);
	CUDA_CHECK();
	CUDA_CALL(cudaMemcpy(
			&ans, a_dev, sizeof(double), cudaMemcpyDeviceToHost));
finally:
	cudaFree(a_dev);
	return ans;
}

一方，ArrayFire で書くとこんな感じ．

https://github.com/hnagata/af-performance/blob/master/af_sum.cxx

#include <vector>
#include <arrayfire.h>
#include "perform.h"

double sum(const std::vector<double> &a) {
	af::array a_dev(a.size(), a.data());
	return af::sum<double>(a_dev);
}

実質1行で済んでしまう．素晴らしい．

パフォーマンス測定には，長さ 1,000,000 の [0, 1] 一様乱数ベクトルを入力として，sum を 1,000 繰り返したトータルタイムを取った．これを3回繰り返して出た最速タイムが次の通り．

	CUDA	ArrayFire
Total time (sec)	3.467	3.042

なんと ArrayFire の方が速い．素人がとりあえず組んだ CUDA よりは，ArrayFire の方がよっぽど最適化されていることがわかる．

合併標本分散

sum はちょっと簡単すぎるので，もう少し面白い例で実験したい．そこで今回は，正規分布の平均差の検定でよく使われる合併標本分散 (pooled sample variance) を取り上げる． 2標本 $X = (X_1, \dots, X_m)$ ， $Y = (Y_1, \dots, Y_n)$ が与えられたとき，合併標本分散 $S^2$ は次の式で計算される．

$\begin{align*} S^2 &= \frac{1}{m + n - 2} \left\{ (m - 1) S_X^2 + (n - 1) S_Y^2 \right\} \\ &= \frac{1}{m + n - 2} \left\{ \sum_{i=1}^m (X_i - \bar X)^2 + \sum_{i=1}^n (Y_i - \bar Y)^2 \right\} \end{align*}$

ただし， $S_X^2, S_Y^2$ は $X, Y$ の標本分散， $\bar X, \bar Y$ は $X, Y$ の標本平均である．

総和の部分は $X$ と $Y$ で独立に計算できるので，CUDA ではストリームを使って2つの総和を非同期に計算することができる．

https://github.com/hnagata/af-performance/blob/master/cuda_pvar.cu

#include <vector>
#include <cuda.h>
#include "perform.h"
#include "cuda_common.h"

const int GRID_SIZE = 448;
const int BLOCK_SIZE = 32;

__global__ void sum_dev(double *a, int n, double *out) {
	__shared__ double work[BLOCK_SIZE];
	const int tidx = threadIdx.x;
	const int bidx = blockIdx.x;
	const int idx = threadIdx.x + blockIdx.x * blockDim.x;
	const int dim = blockDim.x * gridDim.x;

	work[tidx] = 0;

	for (int i = idx; i < n; i += dim) {
		work&#91;tidx&#93; += a&#91;i&#93;;
	}
	__syncthreads();

	int w = BLOCK_SIZE;
	while (w > 2) {
		w /= 2;
		if (tidx < w) work&#91;tidx&#93; = work&#91;tidx + w&#93;;
		__syncthreads();
	}

	if (tidx == 0) {
		out&#91;bidx&#93; = work&#91;0&#93; + work&#91;1&#93;;
	}
}

__global__ void div_dev(double *a, int b, double *out) {
	out&#91;0&#93; = a&#91;0&#93; / b;
}

__global__ void err2_dev(double *a, int n, double *mean, double *out) {
	const int idx = threadIdx.x + blockIdx.x * blockDim.x;
	const int dim = blockDim.x * gridDim.x;
	for (int i = idx; i < n; i += dim) {
		double d = a&#91;i&#93; - mean&#91;0&#93;;
		out&#91;i&#93; = d * d;
	}
}

double sqsum(
	cudaStream_t &stream, 
	const std::vector<double> &a, double *a_dev, int n, double *work)
{
	double out;
	CUDA_CALL(cudaMemcpyAsync(
			a_dev, a.data(), sizeof(double) * n,
			cudaMemcpyHostToDevice, stream));
	sum_dev<<<GRID_SIZE, BLOCK_SIZE, 0, stream>>>(a_dev, n, work);
	sum_dev<<<1, BLOCK_SIZE, 0, stream>>>(work, GRID_SIZE, work);
	div_dev<<<1, 1, 0, stream>>>(work, n, work);
	err2_dev<<<GRID_SIZE, BLOCK_SIZE>>>(a_dev, n, work, a_dev);
	sum_dev<<<GRID_SIZE, BLOCK_SIZE, 0, stream>>>(a_dev, n, work);
	sum_dev<<<1, BLOCK_SIZE, 0, stream>>>(work, GRID_SIZE, work);
	CUDA_CALL(cudaMemcpyAsync(
			&out, work, sizeof(double), cudaMemcpyDeviceToHost, stream));
	CUDA_CALL(cudaStreamSynchronize(stream));
finally:
	return out;
}

double pvar(const std::vector<double> &a, const std::vector<double> &b) {
	const int m = a.size(), n = b.size();
	cudaStream_t a_stream, b_stream;
	double *a_dev = NULL, *b_dev = NULL;
	double *a_work = NULL, *b_work = NULL;
	double a_sqsum, b_sqsum, ans;
	CUDA_CALL(cudaStreamCreate(&a_stream));
	CUDA_CALL(cudaStreamCreate(&b_stream));
	CUDA_CALL(cudaMalloc(
			(void**) &a_dev, sizeof(double) * (m + n + GRID_SIZE * 2)));
	b_dev = a_dev + m;
	a_work = b_dev + n;
	b_work = a_work + GRID_SIZE;

	a_sqsum = sqsum(a_stream, a, a_dev, m, a_work);
	b_sqsum = sqsum(b_stream, b, b_dev, n, b_work);
	ans = (a_sqsum + b_sqsum) / (m + n - 2);
finally:
	cudaFree(a_dev);
	cudaStreamDestroy(a_stream);
	cudaStreamDestroy(b_stream);
	return ans;
}

一方，ArrayFire ではストリームは見えないので，普通にひとつひとつ計算していくことになる．

https://github.com/hnagata/af-performance/blob/master/af_pvar.cxx

#include <vector>
#include <arrayfire.h>
#include "perform.h"

double pvar(const std::vector<double> &a, const std::vector<double> &b) {
	const int m = a.size(), n = b.size();
	af::array a_dev(m, a.data());
	af::array b_dev(n, b.data());
	double a_mean = af::sum<double>(a_dev) / m;
	double b_mean = af::sum<double>(b_dev) / n;
	double a_sqsum = af::sum<double>(af::pow(a_dev - a_mean, 2));
	double b_sqsum = af::sum<double>(af::pow(b_dev - b_mean, 2));
	return (a_sqsum + b_sqsum) / (m + n - 2);
}

ドキュメントには分散を計算する関数 af::var と2乗を計算する関数 af::pow2 があるのだが，実際やってみると “‘***’ is not a member of ‘af'” だったので，これを回避して実装した．プログラムを見ると分かるが，デバイス上の array とホスト側の値で直接演算できたりして，かなり自由度が高い．

で，sum と同様にベクトルの長さ 1,000,000，反復回数 1,000 で実験すると，次のようになった．

	CUDA	ArrayFire
Total time (sec)	7.240	6.762

やっぱり ArrayFire の方が速い．まじですか．

結論

ArrayFire は単にプログラミングを簡単にするだけでなく，パフォーマンス面でも素人が直接 CUDA で組むよりずっと効率よく高速化できることがわかった．ArrayFire は日本ではまだまだ浸透しているとは言えない状況だが，ポテンシャルは非常に高く，またオープンソース化によって誰でも自由に使えるようになったので，これから流行っていくのかもしれない．

補足

実験環境は次の通り．

CPU	Intel Xeon E5-2650
GPU	NVIDIA Tesla C2075
OS	CentOS 6.6
コンパイラ	gcc 4.8.2
CUDA	CUDA 6.0
ArrayFire	@2015-01-30

Posted: 2015-02-05 18:40

最終更新日時: 2015-02-15 14:05

永田晴久

GPGPU ライブラリ “ArrayFire” を生 CUDA と比較する

ArrayFire について

Reduction (sum) のパフォーマンス

合併標本分散

結論

補足