Thrustを使ってCUDAで最近傍法を実装してみた(サンプルコードあり)
この記事では最近傍法(nearest neighbor algorithm)をCUDAで実装する場合のサンプルを紹介しています。Thrustというライブラリを用いることで、GPUコンピューティング、つまりCUDAでのプログラミングが非常に簡単になります。
そこで、Thrustを使ったCUDAのプログラムを紹介するうえで最近傍法の計算部分(距離の計算)を実装してみました。
Thrustを使わない場合のCUDAの最近傍法
注目すべきポイントはdeviceのメモリの管理の部分です。メモリ管理をcudaMallocやcudaMemcpyを使用しなければいけない点です。この処理は直感的ではないため、CUDAについて知らない人が見たときに少しわかりづらいです。また、メモリ確保を怠るとエラーの原因になってしまうため、面倒です。
そこで、Thrustを使うことでメモリ管理が楽になります。
#include<cuda_runtime.h>
#include <device_launch_parameters.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/execution_policy.h>
#define M 5 //点数 m
__global__ void calcNearestNNGpu(double* x1, double* y1, double* x2, double* y2, double* distance, double* dx, double* dy)
{
int threadNumX = blockIdx.x * blockDim.x + threadIdx.x;
if (threadNumX < M){
dx[threadNumX] = x2[threadNumX] - x1[0];
dy[threadNumX] = y2[threadNumX] - y1[0];
distance[threadNumX] = sqrt(pow(dx[threadNumX], 2) + pow(dy[threadNumX], 2));
}
}
int main(){
//hostの変数
double x1[M]; double x2[M]; double y1[M]; double y2[M];
double dx[M]; double dy[M];
double distance[M];
//deviceの変数
double *x1_d; double *x2_d; double *y1_d; double *y2_d;
double *dx_d; double *dy_d;
double *distance_d;
//deviceのメモリの管理
int devSize = sizeof(double)*M;
cudaMalloc((void**)&x1_d, devSize);
cudaMalloc((void**)&x2_d, devSize);
cudaMalloc((void**)&y1_d, devSize);
cudaMalloc((void**)&y2_d, devSize);
cudaMalloc((void**)&dx_d, devSize);
cudaMalloc((void**)&dy_d, devSize);
cudaMalloc((void**)&distance_d, devSize);
cudaMemcpy(x1_d, x1, devSize, cudaMemcpyHostToDevice);
cudaMemcpy(x2_d, x2, devSize, cudaMemcpyHostToDevice);
cudaMemcpy(y1_d, y1, devSize, cudaMemcpyHostToDevice);
cudaMemcpy(y2_d, y2, devSize, cudaMemcpyHostToDevice);
//ホスト上の変数の初期化
for (int initNum = 0; initNum < M; initNum++){
x1[initNum ] = initNum;
x2[initNum ] = initNum;
y1[initNum ] = initNum;
y2[initNum ] = initNum;
std::cout << "x1: " << initNum << "番目: " << x1[initNum] << std::endl;
std::cout << "x2: " << initNum << "番目: " << x2[initNum] << std::endl;
}
const dim3 block(5, 1);
const dim3 grid(1, 1);
//kernelを実行
calcNearestNNGpu << <grid, block >> > (thrust::raw_pointer_cast(&x1_d[0])
, thrust::raw_pointer_cast(&y1_d[0]), thrust::raw_pointer_cast(&x2_d[0])
, thrust::raw_pointer_cast(&y2_d[0]), thrust::raw_pointer_cast(&distance_d[0])
, thrust::raw_pointer_cast(&dx_d[0]), thrust::raw_pointer_cast(&dy_d[0])
);
cudaDeviceSynchronize;
//距離を代入
cudaMemcpy(distance, distance_d, devSize, cudaMemcpyDeviceToHost);
//確認用 距離を表示
for (int j = 0; j < M; j++){
std::cout << "distance: " << j << "番目: " << distance[j] << std::endl;
}
}
Thrustを使う場合のCUDAの最近傍法
Thrustを使う場合、変数の宣言は長くなりますが、Device側(GPU)にHost側(CPU)の値を代入することが簡単になります。 Device側の変数には全て「_d」を付けていますが、x1_dにHost側のx1を代入するだけでDevice側にデータをコピーすることが完了します。
そのため、非常に直感的です。
ただし、1つだけ問題があって、デバッグがしづらいということです。ブレークポイントを設定しても値の参照ができないため、デバッグがしづらいです。
そのため、一々printfなどで値を直接コンソール上に表示したりする必要があります。
#include<cuda_runtime.h>
#include <device_launch_parameters.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/execution_policy.h>
#define M 5 //点数 m
__global__ void calcNearestNNGpu(double* x1, double* y1, double* x2, double* y2, double* distance, double* dx, double* dy)
{
int threadNumX = blockIdx.x * blockDim.x + threadIdx.x;
if (threadNumX < M){
dx[threadNumX] = x2[threadNumX] - x1[0];
dy[threadNumX] = y2[threadNumX] - y1[0];
distance[threadNumX] = sqrt(pow(dx[threadNumX], 2) + pow(dy[threadNumX], 2));
}
}
int main(){
//ホスト上のメモリ定義
thrust::host_vector<double> x1(M);
thrust::host_vector<double> y1(M);
thrust::host_vector<double> x2(M);
thrust::host_vector<double> y2(M);
thrust::host_vector<double> distance(M);
thrust::host_vector<double> dx(M);
thrust::host_vector<double> dy(M);
//デバイス上のメモリ定義
thrust::device_vector<double> x1_d(M);
thrust::device_vector<double> y1_d(M);
thrust::device_vector<double> x2_d(M);
thrust::device_vector<double> y2_d(M);
thrust::device_vector<double> distance_d(M);
thrust::device_vector<double> dx_d(M);
thrust::device_vector<double> dy_d(M);
//host側のデータを1,2,3,4,5で初期化 あらかじめdevice側を初期化してもOK
thrust::sequence(x1.begin(), x1.end());
thrust::sequence(y1.begin(), y1.end());
thrust::sequence(x2.begin(), x2.end());
thrust::sequence(y2.begin(), y2.end());
//ホスト上の変数の初期化
for (int initNum = 0; initNum < M; initNum++){
std::cout << "x1: " << initNum << "番目: " << x1[initNum] << std::endl;
std::cout << "x2: " << initNum << "番目: " << x2[initNum] << std::endl;
}
//デバイスにホストの変数のコピー
x1_d = x1;
y1_d = y1;
x2_d = x2;
y2_d = y2;
dx_d = dx;
dy_d = dy;
const dim3 block(5, 1);
const dim3 grid(1, 1);
//kernelを実行
calcNearestNNGpu << <grid, block >> > (thrust::raw_pointer_cast(&x1_d[0])
, thrust::raw_pointer_cast(&y1_d[0]), thrust::raw_pointer_cast(&x2_d[0])
, thrust::raw_pointer_cast(&y2_d[0]), thrust::raw_pointer_cast(&distance_d[0])
, thrust::raw_pointer_cast(&dx_d[0]), thrust::raw_pointer_cast(&dy_d[0])
);
cudaDeviceSynchronize;
//距離を代入
distance = distance_d;
//確認用 距離を表示
for (int j = 0; j < M; j++){
std::cout << "distance: " << j << "番目: " << distance[j] << std::endl;
}
}
最近傍法におけるThrust部分のまとめ
//ホスト上のメモリ定義
thrust::host_vector<double> x1(M);
//デバイス上のメモリ定義
thrust::device_vector<double> x1_d(M);
//host側のデータを1,2,3,4,5で初期化 あらかじめdevice側を初期化してもOK
thrust::sequence(x1.begin(), x1.end());
//デバイスにホストの変数のコピー
x1_d = x1;
//kernelを実行
calcNearestNNGpu << <grid, block >> > (thrust::raw_pointer_cast(&x1_d[0])
, thrust::raw_pointer_cast(&y1_d[0]), thrust::raw_pointer_cast(&x2_d[0])
, thrust::raw_pointer_cast(&y2_d[0]), thrust::raw_pointer_cast(&distance_d[0])
, thrust::raw_pointer_cast(&dx_d[0]), thrust::raw_pointer_cast(&dy_d[0])
);
//距離を代入
distance = distance_d;
まとめ
Thrustを利用すると、問題点はありますが、非常にプログラミングが楽になります。 そのため、必要に応じて使用してみるといいでしょう。Thrust自体に最小値や最大値を求める機能もあるため、最近傍法との相性も抜群です。 GPUコンピューティングでは最小値を求める処理を自分で書くのは難しいですが、Thrustを使えば簡単です。
今回は距離の計算部分のみでしたが、Thrustを使うことで最小値の計算を1行で終わらせることができます。
