読者です 読者をやめる 読者になる 読者になる

Thrustを使ってCUDAで最近傍法を実装してみた(サンプルコードあり)

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行で終わらせることができます。