cuSOLVER: 鶴亀算を解く

「つる と かめ があわせて３匹います。
　足の数はあわせて10本でした。
　問1: つる と かめ はそれぞれ何匹ですか?
　問2: 足の数があわせて8本なら、それぞれ何匹ですか?」

小学校の算数でやりましたよね、鶴亀算。
コ難しくいえば 連立一次方程式の解を求める問題です。

行列で表現するとこんな。

与えられた A と B から X を求めよ、ってわけ。

CUDA 7.0からこんな問題をCUDAで解いてくれるライブラリ: cuSOLVER がついてます。早速使ってみましたよ。

#include <cuda_runtime.h>
#include <cusolverDn.h> // dense LAPACK

#include <cassert>
#include <iostream>
using namespace std;

template<typename T>
inline size_t bytesof(unsigned int s) { return s * sizeof(T); }

template<typename T>
T* allocate(unsigned int size) {
  T* result = nullptr;
  cudaError_t status = cudaMalloc(&result, bytesof<T>(size));
  assert( status == cudaSuccess );
  return result;
}

int main() {
  cusolverStatus_t status;

  // dense LAPACK
  cusolverDnHandle_t handle;
  status = cusolverDnCreate(&handle);
  assert( status == CUSOLVER_STATUS_SUCCESS );

  int n = 2; // 2x2 正方行列

  float A[] = {  1.0f,  1.0f, 
                 2.0f,  4.0f };
  float* dA = allocate<float>(n*n);
  cudaMemcpy(dA, A, bytesof<float>(n*n), cudaMemcpyHostToDevice);
  int lda = 2;

  // 必要なバッファ量を求め、確保する
  int worksize;
  status = cusolverDnSgetrf_bufferSize(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             &worksize);
  assert( status == CUSOLVER_STATUS_SUCCESS );
#ifdef _DEBUG
  cout << "worksize = " << worksize << endl;
#endif

  float* workspace = allocate<float>(worksize);

  // 計算結果に関する情報
  int* devInfo = allocate<int>(1);
  // ピボット
  int* pivot = allocate<int>(n);

  // LU分解 : dAに結果が求まる(それとpivot)
  status = cusolverDnSgetrf(
             handle,
             n,   // 行
             n,   // 列
             dA,  // A
             lda, // Aのヨコハバ
             workspace,
             pivot,
             devInfo);
#ifdef _DEBUG
  int info;
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  //         鶴と亀の総数, 足の総数
  float B[] = {  3.0f,       10.0f , 
                 3.0f,        8.0f };
  int nrhs = 2; // 問題数
  float* dB = allocate<float>(n*nrhs);
  cudaMemcpy(dB, B, bytesof<float>(n*nrhs), cudaMemcpyHostToDevice);
  int ldb = 2;

  // AX = B を解く (解XはBをoverrideする)
  status = cusolverDnSgetrs(
             handle,
             CUBLAS_OP_T,
             n,     // 行(=列)
             nrhs,  // 問題数
             dA,    // A
             lda,   // Aのヨコハバ
             pivot, // LU分解で得られたピボット
             dB,    // B
             ldb,   // Bのヨコハバ
             devInfo);
#ifdef _DEBUG
  cudaMemcpy(&info, devInfo, sizeof(int), cudaMemcpyDeviceToHost);
  cout << "info = " << info << endl;
#endif
  assert( status == CUSOLVER_STATUS_SUCCESS );

  // 結果を取得し、出力する
  float X[16];
  cudaMemcpy(X, dB, bytesof<float>(n*nrhs), cudaMemcpyDeviceToHost);
  for ( int i = 0; i < nrhs; ++i ) {
    float* q = B + i*2;
    float* a = X + i*2;
    cout << "総数 = " << q[0] << ", 足の数 = " << q[1]
         << "\t 解 :　"
         <<"鶴 = " << a[0] << " , 亀 = " << a[1] <<endl;
  }

  cudaFree(workspace);
  cudaFree(dA);
  cudaFree(dB);
  cudaFree(devInfo);
  cudaFree(pivot);

  cusolverDnDestroy(handle);
}

実行結果:

(original: 2016-06-03 #40)

streamで速くする(3) ～からくり

CPUとGPUはPCI-busに隔てられてそれぞれが勝手に動くことができます。
CPUはGPUの仕事の完了を待って次の仕事を依頼する必要はないんですわ(そうでないとGPUが重たい仕事してる間CPUがぼーっと待ってにゃならんですから)。仕事の完了を待たずに次の仕事を叩き込むことができるってことは、仕事の待ち行列が用意されてるってことで、それがstreamです。

streamに溜まった仕事(メモリ・コピー と kernel実行)はふたつのエンジンが捌きます。kernel実行をCUDA-core群に割り当てるGigaThread Engineと、メモリ・コピーの依頼を受けてPCI-bus経由でCPU/GPU間のデータ転送を司るCopy Engine、この二つは独立して動けます。

んだから、kernel実行中には次のkernel実行は待たされ、メモリ・コピー中は次のメモリ・コピーは待たされるけど、kernel実行とメモリ・コピーとはそれぞれ異なるEngineで処理されるために同時にやれるってわけですわ。

GeForce GTX9xx, GTX10xx になるとCPU→GPU用 と GPU→CPU用のふたつのCopy Engineを持ってます。PCI-busは上り/下りを同時に転送できるのでより速くなります。

GTX980でのtimelineはこんなカンジ:

緑の帯がCPU→GPU、紫の帯がGPU→CPU 両者がoverlapしてる様子が見て取れます。

streamで速くする(2)

streamはその名の通り"一連の(処理の)流れ"です。 streamに乗せられたjob(仕事)は並んだ順に実行されます。 複数のstreamにjobを乗せておけば、GPUはその中から実行可能なものを一つ選んで処理します。 ホスト・メモリ間のデータ転送にGPUは関与しないので、データ転送の間でも他のstreamに乗ったjobを処理できます。

ここまでにいくつかサンプル・コードを紹介してきましたが、それらの中には stream を明示的に使ったものはありません。 けども実は暗黙の default stream が一本だけあって、その default stream に cudaMemcpy やら kernel呼び出しやらのjobを乗っけていたんです。

３本のstream(ssin, scos, stan)を用意し、sin,cos,tanの表作成をそれぞれの stream に乗せて処理しましょう。

#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); }

inline cudaStream_t create_stream() 
  { cudaStream_t s = nullptr; cudaStreamCreate(&s); return s; };

inline void destroy_stream(cudaStream_t stream)
  { cudaStreamDestroy(stream); }

#include <iostream>

using namespace std;

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

void multiple_stream() {
  cout << "--- multiple stream ---" << endl;

  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;
  }

  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;

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

  cudaStream_t ssin = create_stream();
  cudaStream_t scos = create_stream();
  cudaStream_t stan = create_stream();

  cudaEventRecord(start);
  cudaMemcpy(dangle, angle, BN, cudaMemcpyHostToDevice);
  cudaEventRecord(mid);

  kernel_sin<<<grid,block,0,ssin>>>(dsin, dangle, size);
  cudaMemcpyAsync(sin, dsin, BN, cudaMemcpyDeviceToHost, ssin);

  kernel_cos<<<grid,block,0,scos>>>(dcos, dangle, size);
  cudaMemcpyAsync(cos, dcos, BN, cudaMemcpyDeviceToHost, scos);

  kernel_tan<<<grid,block,0,stan>>>(dtan, dangle, size);
  cudaMemcpyAsync(tan, dtan, BN, cudaMemcpyDeviceToHost, stan);

  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 stream : { ssin,  scos, stan }) 
    destroy_stream(stream);
  for ( auto host   : { angle,  sin,  cos,  tan  }) 
    free_host(host);
  for ( auto device : { dangle, dsin, dcos, dtan }) 
    free_device(device);

}

int main() {
  multiple_stream();
}

default stream １本で処理した時のタイムラインがコレ。

対して３本(+default)のstreamでやったのがコチラ。

ね、データ転送とkernel実行が同時に行われてるでしょ。重なった分だけ速くなるんです。

実装上のキモはふたつ:

• cudaMemcpyじゃなく cudaMemcpyAsyncを使うべし。 後者はコピー完了を待たずに(job投入後ただちに)返って来るので、すぐさま次のjobを投入できます。

• ホスト側メモリを malloc/free や new/delete じゃなく cudaMallocHost/cudaFreeHost で確保/解放すべし。 そうしないとデータ転送中に kernelが動いてくれないのですよ。

... streamのおハナシはまだまだ続きます。

(original: 2015-04-16)

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 をお忘れなく。