streamで速くする(1) ～イントロ

三角関数表を作りましょう。 ひと回り360度を百万等分し、360マイクロ度刻みの sin, cos, tan を求めます。

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

// 時間稼ぎ
__device__ float time_waste() {
  float result = 0.0f;
  const int TIMES = 10U;
  for ( int i = 0; i < TIMES; ++i ) {
    result += 1.0f / TIMES;
  }
  return result;
}

// 度→ラジアン
__device__ __forceinline__ float deg2rad(float angle) {
  return angle / 180.0f * 3.01416f * time_waste();
}

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

// cosine
__global__ void kernel_cos(float* out, const float* angle, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = cosf(deg2rad(angle[i]));
  }
}

// tangent
__global__ void kernel_tan(float* out, const float* angle, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = tanf(deg2rad(angle[i]));
  }
}

// helper funcs
template<typename T> inline T* allocate_host(size_t size)
  { T* p = nullptr; cudaMallocHost(&p, size*sizeof(T)); return p; }

template<typename T> inline void free_host(T* host)
  { cudaFreeHost(host); }

template<typename T> inline T* allocate_device(size_t size) 
  { T* p = nullptr; cudaMalloc(&p, size*sizeof(T)); return p; }

template<typename T> inline void free_device(T* device)
  { cudaFree(device); }

inline cudaEvent_t create_event() 
  { cudaEvent_t e = nullptr; cudaEventCreate(&e); return e; };

inline void destroy_event(cudaEvent_t event)
  { cudaEventDestroy(event); }

#include <iostream>

using namespace std;

const unsigned int N = 1000U * 1000U;
const unsigned int BN = N*sizeof(float);

void single_stream() {
  cout << "--- single stream ---" << endl;

  // host memories
  float* angle = allocate_host<float>(N); // in
  float* sin   = allocate_host<float>(N); // out
  float* cos   = allocate_host<float>(N); // out
  float* tan   = allocate_host<float>(N); // out

  for ( unsigned int i = 0; i < N; ++i ) {
    angle[i] = 360.0f * i / N;
  }

  // device memories
  float* dangle = allocate_device<float>(N);
  float* dsin   = allocate_device<float>(N);
  float* dcos   = allocate_device<float>(N);
  float* dtan   = allocate_device<float>(N);

  unsigned int size  = N;
  unsigned int block = 128U;
  unsigned int grid  = (size + block -1U) / block;

  // events
  cudaEvent_t start = create_event();
  cudaEvent_t mid   = create_event();
  cudaEvent_t stop  = create_event();

  cudaEventRecord(start);
  // copy angle[] from host to device
  cudaMemcpy(dangle, angle, BN, cudaMemcpyHostToDevice);
  cudaEventRecord(mid);

  // compute sine
  kernel_sin<<<grid,block>>>(dsin, dangle, size);
  // write-back to host
  cudaMemcpy(sin, dsin, BN, cudaMemcpyDeviceToHost);

  // compute cosine
  kernel_cos<<<grid,block>>>(dcos, dangle, size);
  // write-back to host
  cudaMemcpy(cos, dcos, BN, cudaMemcpyDeviceToHost);

  // compute tangent
  kernel_tan<<<grid,block>>>(dtan, dangle, size);
  // write-back to host
  cudaMemcpy(tan, dtan, BN, cudaMemcpyDeviceToHost);

  cudaEventRecord(stop);
  cudaDeviceSynchronize();

  float elapsed;
  cudaEventElapsedTime(&elapsed, start, mid);
  cout << elapsed << " + ";
  cudaEventElapsedTime(&elapsed, mid, stop);
  cout << elapsed << " = ";
  cudaEventElapsedTime(&elapsed, start, stop);
  cout << elapsed << " [ms]" << endl;
  for ( unsigned int i = 0; i < 5; ++i ) {
    cout << angle[i] << " :"
         << " sin=" << sin[i] 
         << " cos=" << cos[i] 
         << " tan=" << tan[i] 
         << endl; 
  }
  cout << "..." << endl;

  for ( auto event  : { start, mid,  stop }) 
    destroy_event(event);
  for ( auto host   : { angle,  sin,  cos,  tan  }) 
    free_host(host);
  for ( auto device : { dangle, dsin, dcos, dtan }) 
    free_device(device);

}

int main() {
  single_stream();
}

なんともオーソドックスなコードです。入力(angle)をホストからデバイスに送り、"kernel_うんちゃら で計算して、デバイスからホストに書き戻す"を三回やってるだけ。

CUDAを使う限り、ホスト・デバイス間のデータ転送は避けられません。そのことが大きなボトルネックになるのは事実ですが、だからってこれ以上速くはならんかというと、そうでもない。

上記コードでは、ホスト・デバイス間のデータ転送をやってる間、GPUは仕事せずデータ転送が終わるのを待ってます。

ピザが焼け、バイトの兄ちゃんがバイク飛ばして配達に行って、帰って来るまでピザ窯空けたままぼーっと待ってるあほゥなピザ職人はいないでしょ？ 窯が空いたら次の注文受けて焼いとけば、兄ちゃんが配達から帰って来次第すぐさま次の配達できるやないですか。

...てなオハナシを次のエントリで。

(original: 2015-04-15)

Zero Copy

CUDA C Best Practices Guide 9.1.3 に Zero Copy てのが出てきます...なにコレ? ってんで調べてみました。

サンプル: float列の各要素に対し、その平方根をを求める処理を書きました:

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

// out[i] = √in[i] where i = 0..size-1
__global__ void kernel_square_root(const float* in, float* out, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
     out[i]  = sqrtf(in[i]);
  }
}

#include <vector>
#include <numeric>
#include <iostream>

void normal() {
  using namespace std;
  const unsigned int N = 10000U;

  // host領域確保
  float* in  = new float[N];
  float* out = new float[N];

  // device領域確保
  size_t before, after, total;
  cudaMemGetInfo(&before, &total);
  cout << "normal   : device free = " << before << " -> ";
  float*  din; cudaMalloc(&din,  2*N*sizeof(float));
  float* dout; cudaMalloc(&dout, N*sizeof(float));
  cudaMemGetInfo(&after, &total);
  cout << after << " (" << (before-after) << " consumed)" << endl;

  // in[] = { 0.0, 1.0, 2.0, 3.0 ... } 
  iota(in, in+N, 0.0f);

  // host -> device
  cudaMemcpy(din, in, N*sizeof(float), cudaMemcpyHostToDevice);
  // kernel-call
  unsigned int block = 256;
  unsigned int grid = (N + block -1U) / block;
  kernel_square_root<<<grid,block>>>(din, dout, N);
  // device -> host
  cudaMemcpy(out, dout, N*sizeof(float), cudaMemcpyDeviceToHost);

  for ( unsigned int i : { 2, 3, 200, 300} ) {
    cout << out[i] << endl;
  }

  // device領域解放
  cudaFree(din);
  cudaFree(dout);

  // host領域解放
  delete[] in;
  delete[] out;
}

いつものダンドリです:

• ホスト/デバイスの双方に領域を確保
• ホスト領域に入力データを用意
• ホストからデバイスへ(入力データを)コピー
• カーネル呼び出し
• デバイスからホストへ(出力データを)コピー
• ホスト/デバイス領域を解放

一方 zero copy ではデバイス側に領域を確保しません。 そのかわり、デバイスから読み書きできる領域をホスト側に確保します。 zero copyによる「float列の平方根」はこんなコードになります:

void zero_copy() {
  using namespace std;
  const unsigned int N = 10000U;

  // host領域確保
  size_t before, after, total;
  cudaMemGetInfo(&before, &total);
  cout << "zero-copy: device free = " << before << " -> ";
  float* in;  cudaHostAlloc(&in,  N*sizeof(float), cudaHostAllocMapped);
  float* out; cudaHostAlloc(&out, N*sizeof(float), cudaHostAllocMapped);

  // hostにマップされたdeviceポインタを取得
  float* din;  cudaHostGetDevicePointer(&din,  (void*)in,  0);
  float* dout; cudaHostGetDevicePointer(&dout, (void*)out, 0);

  cudaMemGetInfo(&after, &total);
  cout << after << " (" << (before-after) << " consumed)" << endl;

  // in[] = { 0.0, 1.0, 2.0, 3.0 ... } 
  iota(in, in+N, 0.0f);

  // kernel-call
  unsigned int block = 256;
  unsigned int grid = (N + block -1U) / block;
  kernel_square_root<<<grid,block>>>(din, dout, N);
  // kernel処理完了を待つ
  cudaStreamSynchronize(nullptr);

  for ( unsigned int i : { 2, 3, 200, 300} ) {
    cout << out[i] << endl;
  }

  // host領域解放
  cudaFreeHost(in);
  cudaFreeHost(out);
}

違いは３つ:

• ホスト側の領域確保には cudaHostAlloc( ... cudaHostAllocMapped) を用いる。
• デバイス側の領域確保は行わず、cudaHostGetDevicePointer でホスト側領域に対応するデバイス・ポインタを手に入れ、kernelに渡す。
• kernel呼び出しのあと、kernelの処理が完了するのを待つ(cudaMemcpyは処理完了を暗黙裡に待つが、それが行われないため)。

両者を実行しました:

int main() {
  // zero copy できるか確認
  cudaDeviceProp prop;
  cudaGetDeviceProperties(&prop, 0);
  if ( !prop.canMapHostMemory ) {
    std::cerr << "sorry, not supported." << std::endl;
    return 0;
  }
  // zero copy を有効にする
  cudaSetDeviceFlags(cudaDeviceMapHost);

  normal();
  zero_copy();
}

ご覧のとおり、zero copy版はデバイス・メモリを消費していません。

cudaHostAlloc( ... cudaHostAllocMapped) で確保されたホスト側領域は物理アドレスが固定(ピン留め)されていて、kernelはこの領域に対し、PCI-bus越しに読み書きするってからくりです。

読み書きが発生するたびにPCI-busを行き来しますから、kernelがこの領域に対して何度も読み書きするとか、(連続領域をリニアに、ではなく)飛び飛びのアドレスにアクセスすると、kernel-callをcudaMemcpyで挟むいつもの手順に比べてパフォーマンスを著しく損なうことになります。

Tegra K1/X1 なんてな integrated-GPU だとCPUとGPUはひとつのメモリ空間を共有してるんで cudaMemcpyは無駄であり、ZeroCopyが効果的に機能しますですよきっと。

(original: 2015-08-03)

後付けピン留め

cudaMallocHostを使ってピン留めされたHost-memoryを確保すればcudaMemcpyのスキマを潰すことができ、その分速くなる(てかそうしないと複数streamでのoverlapができん)のですが、あらかじめ(mallocやnewで)確保された領域を後付けでピン留め...できるんですわこれが。

cudaHostRegisterってAPIがありまして、領域の先頭アドレスと大きさ(バイト数)を与えることでピン留めしてくれます。ピンを外すのはcudaHostUnregisterで。

#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 HOSTREGISTER

int main() {

  const unsigned int N = 100000;

// Host-memoryを確保
  float* h_angle;
  float* h_sine;
  h_angle = (float*)malloc(N*sizeof(float));
  h_sine  = (float*)malloc(N*sizeof(float));
#if defined(HOSTREGISTER)
　　// 後付けピン留め
  cudaHostRegister(h_angle, N*sizeof(float), cudaHostRegisterDefault);
  cudaHostRegister(h_sine , N*sizeof(float), cudaHostRegisterDefault);
#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(HOSTREGISTER)
  // ピン留め解除
  cudaHostUnregister(h_angle);
  cudaHostUnregister(h_sine );
#endif
  free(h_angle);
  free(h_sine );
  cudaDeviceReset();
}

timelineはこんな↓、スキマ詰まってますね。

ただね、このcudaHostRegister使った後付けピン留めってかなり時間がかかるみたい。なので頻繁に留め/外しを繰り返すもんじゃなさそうです。

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

さっきのarticle、総計を求めるためだけに thrust::reduce を使いましたけど、Thrustはすっごく便利で使いでがあるんですよ。

たとえば thrust::device_vector<T> は cudaMalloc/cudaFree で取得/解放するデバイス・メモリーが可変長配列であるかのようにふるまってくれるクラス。 デバイス・メモリーがstd::vector<T> みたいに使えちゃう。

「モンテカルロで円周率」を Thrust色濃いメにre-writeしてみました。
ぱっと見これがCUDAとは思えぬコードができましたよ♪

/*
  DO NOT FORGET OPTION: --expt-extended-lambda 
 */
// CUDA Runtime
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

// cuRAND
#include <curand.h>

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

#include <iostream>
using namespace std;

int main() {

  const size_t N = 1000 * 1000;

  thrust::device_vector<float> x(N);
  thrust::device_vector<float> y(N);
  thrust::device_vector<float> count(N);

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

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

  // i番目の点(x[i],y[i])が1/4円の中に入っていれば count[i] = 1
  cout << "--- is the point in the circle? ---" << endl;
  thrust::transform(begin(x), end(x), begin(y), begin(count), 
                    [] __device__ (float x, float y) { 
                      return ( x*x + y*y < 1.0f ) ? 1.0f : 0.0f; 
                    }); 

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

  curandDestroyGenerator(gen);

}

(original: 2015-04-06)
追記: CUDA 7.5 以降、device-lambda が使えるようになりました。
コンパイル・オプション: --expt-extended-lambda をお忘れなく。

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

二次元平面上、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)