出典: トータル・ディスクロージャ・サイト（事実をありのままに）

モンテカルロ法などに効果を発揮する、CUDA Toolkit付属の乱数ライブラリ、CURAND。CUDA 4.1では、このCURANDに、MRG32k3aとMersenne Twisterの、2種類のアルゴリズムが新たに加わった。そこで、今回は、これらの新アルゴリズムを試してみた。

目次 1 テストの内容 1.1 CUDA Cプログラム 1.2 実験環境 1.3 実行結果 2 まとめ

テストの内容

モンテカルロ法を用いて、円周率の値を求めるプログラムを、CURANDの各種アルゴリズムで実装し、実行時間と精度の面から比較を行う。

CUDA Cプログラム

#include <cstdio>

#include <sys/time.h>

#include <curand.h>
#include <thrust/device_ptr.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/reduce.h>

struct func : thrust::unary_function<int, int> {
  double *d_a;
  func(double *d_a_) : d_a(d_a_)
  {}
  __device__ int
  operator()(int i) const
  {
    double x = d_a[i * 2 + 0];
    double y = d_a[i * 2 + 1];
    double r2 = x * x + y * y;
    return (r2 <= 1.0) ? 1 : 0;
  }
};

void
benchmark(const size_t &n, double *d_a, const curandRngType_t &rngtype,
          const char *name)
{
  struct timeval tv0, tv1;
  gettimeofday(&tv0, NULL);
  curandGenerator_t generator;
  curandCreateGenerator(&generator, rngtype);
  curandGenerateUniformDouble(generator, d_a, n * 2);
  thrust::counting_iterator<int> count(0);
  thrust::transform_iterator<func, thrust::counting_iterator<int> >
    trans(count, func(d_a));
  int sum = thrust::reduce(trans, trans + n);
  double pi = 4.0 * double(sum) / double(n);
  curandDestroyGenerator(generator);
  gettimeofday(&tv1, NULL);
  double error = std::abs(pi - M_PI) / M_PI;
  int sec = tv1.tv_sec - tv0.tv_sec;
  int usec = tv1.tv_usec - tv0.tv_usec;
  if (usec < 0) {
    --sec;
    usec += 1000000;
  }
  printf("%s (%d Points)\nPI:\t%.70f\nError:\t%.70f\nTime:\t%d.%06d\n",
         name, n, pi, error, sec, usec);
}

int
main()
{
  const size_t n = 300000000;

  double *d_a;
  cudaMalloc((void **)&d_a, n * 2 * sizeof(double));
  benchmark(n, d_a, CURAND_RNG_PSEUDO_XORWOW, "XORWOW");
  benchmark(n, d_a, CURAND_RNG_PSEUDO_MRG32K3A, "MRG32k3a");
  benchmark(n, d_a, CURAND_RNG_PSEUDO_MTGP32, "Mersenne Twister");
  benchmark(n, d_a, CURAND_RNG_QUASI_SOBOL64, "Sobol64");
  benchmark(n, d_a, CURAND_RNG_QUASI_SCRAMBLED_SOBOL64, "Scrambled Sobol64");
  cudaFree(d_a);
  return 0;
}

Thrustテンプレートを使っているため、見慣れないテンプレート関数が出てきているが、XとY座標に、それぞれ0 - 1の倍精度一様乱数を3億要素生成し、二乗和が1以下の場合に1、1を超えた場合に0を返すフィルタを通し、その結果の和を取って、4を掛けて要素数3億で割るプログラムである。後述するように、実行環境がメモリ6 GBのTesla M2090であるので、要素数を3億に設定しているが、メモリの少ないGPUによる実行の場合は、要素数nを少なくして計算することが望ましい。

実験環境

CPU:	Six-Core Xeon X5690 3.46 GHz × 2
GPU:	Tesla M2090 × 4
OS:	CentOS 6.0
CUDA:	CUDA 4.1 RC1

実行結果

XORWOW (300000000 Points)
PI:	3.1415300266666665862658192054368555545806884765625000000000000000000000
Error:	0.0000199347687724466886792767245095703287915966939181089401245117187500
Time:	0.110700
MRG32k3a (300000000 Points)
PI:	3.1416060933333334048711549257859587669372558593750000000000000000000000
Error:	0.0000042780032366486873097945371757777621724017080850899219512939453125
Time:	0.255586
Mersenne Twister (300000000 Points)
PI:	3.1414099733333333830387346097268164157867431640625000000000000000000000
Error:	0.0000581489316417232939178394435231211900827474892139434814453125000000
Time:	0.207856
Sobol64 (300000000 Points)
PI:	3.2915026266666664866988867288455367088317871093750000000000000000000000
Error:	0.0477178264679146910132168102336436277255415916442871093750000000000000
Time:	0.152308
Scrambled Sobol64 (300000000 Points)
PI:	3.9815186533333335461293245316483080387115478515625000000000000000000000
Error:	0.2673567493811729312014335846470203250646591186523437500000000000000000
Time:	0.152550

MRG32k3aは、従来手法に比べて、実行時間はかなり長いが、精度はかなり向上している。Mersenne Twisterに関しては、今回は実行時間でも精度でも、従来手法のXORWOWに負けているが、アプリケーションによっては有利な局面も出てくることも考えられる。

まとめ

CUDA 4.1に実装された、CURANDの新アルゴリズムを使うことによって、より精度の高いモンテカルロ法が、実現可能になると考えられる。