モンテカルロ法で円周率を求める

二次元平面上、0 < x <= 1, 0 < y <= 1 の領域に、デタラメに点を打ちます。その点と原点(0,0)との距離が1以下であるなら、点は1/4円の内部にあります。

大量の点　(x[i],y[i]) : i = 0, 1...N-1 をデタラメに打ち、1/4円の中に入った点の数をPとすると、P/N はだいたい π/4 になるハズね。CUDAでやってみましょうね。

// CUDA Runtime
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

// cuRAND
#include <curand.h>

// Thrust
#include <thrust/device_ptr.h>
#include <thrust/reduce.h>

#include <iostream>
using namespace std;

__global__ void kernel_is_in_circle(const float* x, const float* y, 
                                    float* count, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    count[i] = ( x[i]*x[i] + y[i]*y[i] < 1.0f ) ? 1.0f : 0.0f; 
  }
}

int main() {
  const unsigned int N = 1000 * 1000;

  float* x = nullptr;
  float* y = nullptr;
  float* count = nullptr;

  cudaMalloc(&x,     N*sizeof(float));
  cudaMalloc(&y,     N*sizeof(float));
  cudaMalloc(&count, N*sizeof(float));

  curandGenerator_t gen;
  // メルセンヌ・ツイスタ
  curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_MTGP32);
  curandSetPseudoRandomGeneratorSeed(gen, 1234ULL); // seedはテキトーな値

  // 0～1の一様乱数列を生成する
  cout << "--- generating uniform random numbers ---" << endl;
  curandGenerateUniform(gen, x, N);
  curandGenerateUniform(gen, y, N);

  // i番目の点(x[i],y[i])が1/4円の中に入っていれば count[i] = 1
  cout << "--- is the point in the circle? ---" < endl;
  unsigned int block = 256;
  unsigned int grid = (N + block -1) / block;
  kernel_is_in_circle<<<grid,block>>>(x, y, count, N);

  // count[N]の総計をNで割って4倍すれば π になるハズ
  cout << "--- accumulating the result(pi) ---" << endl;
  thrust::device_ptr<float> dev(count);
  float pi = thrust::reduce(dev, dev+N) * 4.0f / N;
  cout << pi << endl;

  curandDestroyGenerator(gen);

  cudaFree(x);
  cudaFree(y);
  cudaFree(count);

}

百万個の点を打って勘定したら 3.14116、まぁいいセンいってんじゃないでしょか。

(original: 2015-04-06)

CUDAを遊ぶための最小構成

CUDAは「GPUを汎用計算に使っちゃおう」を目的としたNVIDIAのGPUコンピューティングに付けられた名前。当然ながらNVIDIA製GPUを積んだビデオ・カードが必要です。

NVIDIAのページ:CUDA GPUsを眺めるとCUDAが動くGPUが列挙されています。

先日アキバで最安のビデオ・カードを買ってきました、GeForce GT720を積んでて税込6,500円也。グラフィクス性能はお値段から推して知るべしなのですが、ホビーやお試しでCUDAを遊ぶつもりならこれで十分です。(おサイフに余裕があるなら GTX750, 750Ti あたりがオススメです。)

開発環境はWindowsの場合 CUDA Toolkit と Visual Studio、
CUDA Toolkit は NVIDIA から無償でダウンロードできます。

CUDA 最新版 7.0 のRelease noteに、こんな一覧が見つかりました。

無償の Visual Studio 2013 Community が使えます。 Visual Studio のインストールの後、CUDA をインストールしてください。
32bitはダメよとあります。が、どっちみち CUDA 7.0 以降 32bitサポートは限定的になっていて、Toolkitについてくるいろんなライブラリが32bit版では提供されていないのですよ。

...てなわけで、64bitのWindows7/8.1 が動いているなら、あとはCUDA 7.0 と Visual Studio 2013 Community、それとお安いビデオ・カードがあればCUDAで遊べます。

(original: 2015-04-03) 追記:
* 今ならGTX9xx, GTX10xx でしょうけど、それより古いものでも大丈夫。でもあんまり古いのはちょっと... CUDA Toolkit 8.0 以降、Fermi-architectureのGPUはdeprecated扱いになってるみたい。ここで紹介した GT720 はKepler, GTX750(Ti) および GTX9xx はMaxwell、GTX10xx はPascalです。
* CUDA Toolkit 8.0 で Visual Studio 2015 Community に対応しました。

無害な warning を黙らせる

Visual Studio に CUDA Toolkit をインストールし、プロジェクト・テンプレートに追加された CUDA Runtime でひな型プロジェクトを吐かせると、プロジェクト内には小さなコード:kernel.cuがひとつぽつんと作られます。

ともかくもビルドしてみると...exeができてはくれるものの、kernel.cuのコンパイル中に warning: C4819 がざくざく出てきます。

コレ、CUDA Toolkit が提供してるヘッダが暗黙裡に#includeされ、しかもそのヘッダがUTF-8 encodingであることに起因します。

Visual C++ はdefault-encodingを CP932(shift-jis) としてるところに UTF-8 文字が混じるもんだから encoding の不整合を見つけて warning を吐いてるんですな。実害はないとはいえコンパイルのたんびに山ほどの warning を見せられるのは精神衛生上よろしくありません。

この warning を黙らせるには プロジェクト・プロパティ: CUDA C/C++ : Command Line に -Xcompiler -wd4819 を設定します。

(original: 2015-04-03)

ピン留め(Page-locked)されたHost-memoryの効果

CUDAによる一連の処理は、

• cudaMemcpy でデータ(input)をHostからDeviceにコピー
• kernelを実行: inputを読んで処理してoutputに書く
• cudaMemcpy でデータ(output)をDeviceからHostにコピー

ってゆー一連のダンドリになります。
このとき Host-memoryをmalloc/new で確保するのとcudaHostAlloc/cudaMallocHostで確保するのとでは動きがかなり異なります。

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <math_functions.h>

__global__ void kernel_sine(const float* angle, float* sine, unsigned int size) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    sine[i] = sinf(angle[i]);
  }
}

#include <cstdlib>
#include <iostream>

// どちらかひとつを有効にしておくれ
#define MALLOC
// #define CUDAMALLOCHOST

int main() {

  const unsigned int N = 100000;

// Host-memoryを確保
  float* h_angle;
  float* h_sine;
#if defined(MALLOC)
  h_angle = (float*)malloc(N*sizeof(float));
  h_sine  = (float*)malloc(N*sizeof(float));
#elif defined(CUDAMALLOCHOST)
  cudaMallocHost(&h_angle, N*sizeof(float));
  cudaMallocHost(&h_sine , N*sizeof(float));
#endif

  // Device-memoryを確保
  float* d_angle;
  float* d_sine;
  cudaMalloc(&d_angle, N*sizeof(float));
  cudaMalloc(&d_sine , N*sizeof(float));
  
  // HtoD, kernel, DtoH を二回実行
  cudaMemcpyAsync(d_angle, h_angle, N*sizeof(float), cudaMemcpyDefault);
  kernel_sine<<<(N+255)/256,256>>>(d_angle, d_sine, N);
  cudaMemcpyAsync(h_sine , d_sine , N*sizeof(float), cudaMemcpyDefault);

  cudaMemcpyAsync(d_angle, h_angle, N*sizeof(float), cudaMemcpyDefault);
  kernel_sine<<<(N+255)/256,256>>>(d_angle, d_sine, N);
  cudaMemcpyAsync(h_sine , d_sine , N*sizeof(float), cudaMemcpyDefault);

  cudaStreamSynchronize(nullptr); // defaul-stream上の処理が完了するまで待機
  std::cout << "done." << std::endl;

  // あとしまつ
  cudaFree(d_angle);
  cudaFree(d_sine );

#if defined(MALLOC)
  free(h_angle);
  free(h_sine );
#elif defined(CUDAMALLOCHOST)
  cudaFree(h_angle);
  cudaFree(h_sine );
#endif
  cudaDeviceReset();
}

実行時のtimelineはそれぞれこんな。

maloc

cudaMallocHost

mallocだと一回目と二回目との間にスキマができてますが、cudaMallocHostではきっちり詰まってますよね。つまりそれだけ速いってことです。

PCI-bsを介してHost-Device間のコピーが行われるとき、DMA(Direct Memory Access)が使われます。CPUの動作とは無関係に(独立して)勝手にメモリをアクセスする機能です。DMAが機能するためにはコピー元/先の物理アドレスが固定されていなくてはなりません。

malloc/newで確保された領域は仮想アドレス空間にあり、物理アドレスが定まっていないのでDMAが使えない、そこで一旦 Staging-bufferと呼ばれる固定領域にコピーされた後、DMAが動きます。スキマが空くのはこのため。

一方cudaMallocHostはハナっから固定された(Page-locked)領域を確保するのでいきなりDMAが動けるってスンポーです。

複数のStreamを使ったmemory-copyとkernelとのoverlapはHost-memoryが固定されていないとその効果が期待できません。

※ ただし、固定域をたくさん確保するとマシン全体のパフォーマンスを落とすことになりかねません。ご利用は計画的に。

cudaMemcpyのコピー方向とUVA

cudaMemcpy の引数は コピー先, コピー元, バイト数, そして コピー方向。 コピー方向は コピー先/元がHostかDeviceのそれぞれに応じて4種類...ともうひとつ cudaMemcpyDefault てのがあります。

cudaMemcpyDefault はコピー元/先に与えたポインタがHost/Deviceのどちらかを判別し善きに計らってくれます。 が cudaMemcpyDefault を使うには「Unified Virtual Addressing をサポートしていること」て但し書きがあります。

UVA(Unified Virtual Addressing) とは、CPUとGPUのメモリ空間を仮想的にひとつの空間に配置できるてゆーよくわかんない機能でして、これが有効じゃないと与えられたポインタがHost/Deviceのどっちか判別できんとのこと。

で、UVAが有効か否かはこんなコードで確認できます:

#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include <iostream>
using namespace std;

int main() {

    {
    int device;
    cudaGetDevice(&device);
    cudaDeviceProp property;
    // 現デバイスのプロパティを取得し、
    cudaGetDeviceProperties(&property, device);
    // unifiedAddressing != 0 なら UVA有効
    if ( property.unifiedAddressing ) {
      cout << "UVA enabled, cudaMemcpyDefault can be used." << endl;
    } else {
      cout << "sorry, no-UVA" << endl;
    }
  }

  // 試しに cudaMemcpuDefault を使ってみる
  int* host_ptr;   host_ptr = new int[1];
  int* device_ptr; cudaMalloc(&device_ptr, sizeof(int));

  // host->device
  *host_ptr = 12345;
  cudaMemcpy(device_ptr, host_ptr, sizeof(int), cudaMemcpyDefault);
  // device->host
  *host_ptr = 0;
  cudaMemcpy(host_ptr, device_ptr, sizeof(int), cudaMemcpyDefault);

  if ( *host_ptr == 12345 ) {
    cout << "ok." << endl;
  } else {
    cout << "oops!." << endl;
  }

  delete[] host_ptr;
  cudaFree(device_ptr);
  cudaDeviceReset();
}

UVAは近頃の大抵のGPUでサポートしてるのですが、32bitだとダメです。32bitではGPU/CPU双方の空間をひとつにするに十分な大きさじゃないんですな。

東方算程譚：只今熱烈引越中

ココからお引越しですー