まえがき

Unified Memory のおかげで CUDA もだいぶ高級言語らしい書き方ができるようになってきました。ただ、いまいち使われてないようなので典型的な書き方を示しながら、CUDA をおさらいしてみます。

注意点

以下のコードではソースを直感的にわかりやすくするために意識的にエラー処理を省いています。
CUDA でプログラムする局面では double と float、時には int の全てに対応する必要がある場合が多いです。そのためにテンプレートを使っています。

1次元

配列の 2 乗をインプレースでやる例を考えてみます。
・ Unified Memory を使用するには、cudaMallocManaged を呼びます。

最もシンプルな例：

1.cu

# include	<iostream>
using namespace std;

template < typename F > F*
MallocManaged( size_t N ) {
	F*	_;
	cudaMallocManaged( &_, N * sizeof( F ) );
	return _;
}

//	vv	#1
template	< typename F >	__global__	void
SQUARE( F* _ ) {
	_[ threadIdx.x ] *= _[ threadIdx.x ];
}
//	^^

template	< typename F >	void
Main( size_t W ) {

	auto _ = MallocManaged< F >( W );
	for ( size_t x = 0; x < W; x++ ) _[ x ] = x;

//	vv	#2
	SQUARE<<< 1, W >>>( _ );
//	^^
	cudaDeviceSynchronize();
	for ( size_t x = 0; x < W; x++ ) cout << ' ' << _[ x ];
	cout << endl;
	cudaFree( _ );
}
int
main( int argc, char* argv[] ) {
	Main< float >( 4 );
}

実行結果

0 1 4 9

このコードは W がある一定の数値（多くは 1024、以後 maxThreadsPerBlock と呼びます）以下の場合にはうまく動きます。しかし W が maxThreadsPerBlock より多くなると実行時にエラーとなります。
maxThreadsPerBlock は以下のようにしてとることができます。

static	inline	size_t
maxThreadsPerBlock( int ID = 0 ) {
	cudaDeviceProp _;
	cudaGetDeviceProperties( &_, ID );
	return (size_t)_.maxThreadsPerBlock;
}

W が maxThreadsPerBlock を超えても動くようにするために、多くの場合は以下のようにソースを変更します。

1.mod.cu

//	#1
template	< typename F >	__global__	void
SQUARE( F* _, size_t W ) {
	auto x = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
	if ( x >= W ) return;
	_[ x ] *= _[ x ];
}
//	^^
//	#2
	static	const	size_t	BLOCK_SIZE = 256;
	SQUARE<<< ( W + BLOCK_SIZE ) / BLOCK_SIZE, BLOCK_SIZE >>>( _, W );
//	^^

BLOCK_SIZE は maxThreadsPerBlock 以下である必要があります。実際にいくつにするかは経験的な判断になりますが、多くの場合ワープサイズである 32 の倍数です。

このやりかただと、配列の要素数以上のスレッドが走る可能性があるため、不要なスレッドは何もしないようにするために
if ( x >= W ) return;
があります。

struct をカーネルの引数として渡す。

取得した Unified Memory とそのサイズを struct として持っておきたいというような場合、ありますよね。スコープから出た時に cudaFree を必ずよんでくれるようにできますし。そのような場合は以下のようにするといいと思います。

2.cu

# include	<iostream>
using namespace std;

template < typename F > F*
MallocManaged( size_t N ) {
    F*  _;
    cudaMallocManaged( &_, N * sizeof( F ) );
    return _;
}


template	< typename F >	struct
vArray {
	F*		_;
	size_t	w;

	vArray( F* _, size_t w )
	:	_( _ )
	,	w( w ) {
	}
};

template	< typename F >	struct
Array	: vArray< F > {

	~
	Array() {
		cudaFree( vArray< F >::_ );
	}
	Array( size_t w )
	:	vArray< F >( MallocManaged< F >( w ), w ) {
	}
};

template	< typename F >	__global__	void
SQUARE( vArray< F > _ ) {
	auto x = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
	if ( x >= _.w ) return;
	_._[ x ] *= _._[ x ];
}

template	< typename F >	void
Main( size_t W ) {

	Array< F > _( W );
	for ( size_t x = 0; x < W; x++ ) _._[ x ] = x;
	static	const	size_t	BLOCK_SIZE = 256;
	SQUARE<<< ( W + BLOCK_SIZE ) / BLOCK_SIZE, BLOCK_SIZE >>>( _ );
	cudaDeviceSynchronize();
	for ( size_t x = 0; x < W; x++ ) cout << ' ' << _._[ x ];
	cout << endl;
}
int
main( int argc, char* argv[] ) {
	Main< double >( 4 );
}

ところで、何故以下のように単純にしないで、vArray と Array の 2 段構えなんだろう、という疑問がわきませんか？

NG.cu

template	< typename F >	struct
Array {
	F*		_;
	size_t	w;

	~
	Array() {
		cudaFree( vArray< F >::_ );
	}
	Array( size_t w )
	:	_( MallocManaged< F >( w ) )
	,	w( w ) {
	}
};

template	< typename F >	__global__	void
SQUARE( Array< F > _ ) {
	auto x = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
	if ( x >= _.w ) return;
	_._[ x ] *= _._[ x ];
}

実はこうすると、cudaFree が同じアドレスに対して 2 回呼ばれることになります。カーネルに渡した時点で Array はコピーされてしまうからです。これに対処する方法はいくつか考えられますが、この 2 段構えの方法が一番良さそうです。他の方法は以下のように考えられます。

コピーコンストラクタを書いて中でメモリをコピーする。

これはありかもしれないですが、パフォーマンスに影響を与えると思われます。

ポインタ渡しか参照渡しにする

カーネルはポインタで渡された引数をデバイス側のメモリだと思ってしまうのでうまくいきません。参照渡しも同様です。

二次元

二次元の例も示しておきます。

3.cu

# include	<iostream>
using namespace std;

template < typename F > F*
MallocManaged( size_t N ) {
    F*  _;
    cudaMallocManaged( &_, N * sizeof( F ) );
    return _;
}

template	< typename F >	struct
vMatrix {
	F*		_;
	size_t	h;
	size_t	w;

	vMatrix( F* _, size_t h, size_t w )
	:	_( _ )
	,	h( h )
	,	w( w ) {
	}
};

template	< typename F >	struct
Matrix	: vMatrix< F > {

	~
	Matrix() {
		cudaFree( vMatrix< F >::_ );
	}
	Matrix( size_t h, size_t w )
	:	vMatrix< F >( MallocManaged< F >( h * w ), h, w ) {
	}
};

template	< typename F >	__global__	void
SQUARE( vMatrix< F > _ ) {
	auto x = (size_t)blockIdx.x * blockDim.x + threadIdx.x;
	if ( x >= _.w ) return;
	auto y = (size_t)blockIdx.y * blockDim.y + threadIdx.y;
	if ( y >= _.h ) return;
	_._[ y * _.w + x ] *= _._[ y * _.w + x ];
}

template	< typename F >	void
Main( size_t H, size_t W ) {

	Matrix< F > _( H, W );
	for ( size_t y = 0; y < H; y++ ) {
		for ( size_t x = 0; x < W; x++ ) {
			_._[ y * W + x ] = y * W + x;
		}
	}
	static	const	size_t	B_H = 32;
	static	const	size_t	B_W = 32;
	SQUARE<<<
		dim3( ( W + B_W - 1 ) / B_W, ( H + B_H - 1 ) / B_H )
	,	dim3( B_W, B_H )
	>>>( _ );
	cudaDeviceSynchronize();
	for ( size_t y = 0; y < H; y++ ) {
		for ( size_t x = 0; x < W; x++ ) {
			cout << ' ' << _._[ y * W + x ];
		}
		cout << endl;
	}
}
int
main( int argc, char* argv[] ) {
	Main< double >( 2, 3 );
}

実行結果

0 1 4
9 16 25

2次元の場合のブロックサイズは maxThreadsPerBlock が 1024 の場合 32 x 32 の一択な感じです。

最後に

Unified Memory はホストメモリにアクセスするときにどうしてもオーバーヘッドがあるので、例えばホストメモリ上で全要素にシークエンシャルにアクセスしたりすると、えらくパフォーマンスが低下します。なのでそういうアプローチをさけるか Unified Memory を諦めるかの判断が必要な場合が生じるので、パフォーマンスのモニターが必須です。