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

東方算程譚

επιστημηがヨタをこく、弾幕とは無縁のCUDAなタワゴト

on-memory database の検索 : SELECT JOIN ~

SELECT WHERE の応用(?)として SELECT JOIN の試作。

SELECT * 
  FROM A, B
  JOIN A.id = B.id
  WHERE A.weight < 0.1;

ふたつの配列 A[Na], B[Nb] があって、A と B の直積すなわち A[x]. B[y] (x = 0..Na-1, y = 0..Nb-1) のすべてのペア(Na*Nb通り)に対し、A[x].id == B[y].id かつ A[x].weight < 0.1 を満たす x,y の組を列挙する、と。

直積ですから x,y に対する 二重のfor-loop で全組み合わせが作れます。CUDAではスレッド数を決定する grid, block が三次元(dim3)なので三重のloopまでならするっと書けます。

/*
   SELECT * FROM arecord, brecord
     JOIN arecord.id = brecord.id
     WHERE arecord.weight < 0.1
*/

// nvcc fake_DB.cpp --expt-extended-lambda 
// CUDA 8.0

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

struct Arecord {
  int   id;
  float weight;
};

struct Brecord {
  int   id;
  float height;
};

template<typename Join, typename Where>
__global__ void kernel_select(
    const Arecord* arecords, unsigned int asize,
    const Brecord* brecords, unsigned int bsize, 
    Join jpredicate, Where wpredicate,
    int* count, uint2* indices) {
  unsigned int sx = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int sy = blockDim.y * blockIdx.y + threadIdx.y;
  if ( sx ==0 && sy == 0 ) *count = 0;
  __syncthreads();

  for ( unsigned int x = sx; x < asize; x += gridDim.x * blockDim.x ) {
    for ( unsigned int y = sy; y < bsize; y += gridDim.y * blockDim.y ) {
      if ( jpredicate(arecords[x],brecords[y]) &&
           wpredicate(arecords[x],brecords[y]) ) {
        indices[atomicAdd(count, 1)] = uint2{x,y};
      }
    }
    
  }
}

#include <random>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;

int main() {
  const int Na = 100;
  const int Nb = 100;
  Arecord arecords[Na];
  Brecord brecords[Nb];

  // fill records[] with random numbers
  mt19937 gen;
  uniform_int_distribution<int> idist(1,10);
  uniform_real_distribution<float> fdist;
  generate_n(arecords, Na, [&]() { return Arecord{ idist(gen),fdist(gen) }; });
  generate_n(brecords, Nb, [&]() { return Brecord{ idist(gen),fdist(gen) }; });

  // result
  uint2 indices[Na*Nb];
  int   count;

  // allocate device-mem.
  Arecord* dev_arecords;
  cudaMalloc(&dev_arecords, Na*sizeof(Arecord));
  Brecord* dev_brecords;
  cudaMalloc(&dev_brecords, Nb*sizeof(Brecord));
  uint2* dev_indices;
  cudaMalloc(&dev_indices, Na*Nb * sizeof(uint2));
  int* dev_count;
  cudaMalloc(&dev_count, sizeof(int));

  // copy records from HOST to DEVICE
  cudaMemcpy(dev_arecords, arecords, Na*sizeof(Arecord), cudaMemcpyHostToDevice);
  cudaMemcpy(dev_brecords, brecords, Nb*sizeof(Brecord), cudaMemcpyHostToDevice);

  auto join  = [] __device__(const Arecord& arec, const Brecord& brec) { return arec.id == brec.id; };
  auto where = [] __device__(const Arecord& arec, const Brecord& brec) { return arec.weight < 0.1f; };

  dim3 grid { 2, 2 };
  dim3 block { 38, 8 };
  kernel_select<<<grid, block>>>(
    dev_arecords, Na, // FROM arecords
    dev_brecords, Nb, //     ,brecords
    join,             // JOIN arecords.id = brecords.id
    where,            // WHERE arecords.weight < 0.1
    dev_count, dev_indices); // SELECTed int-pair

  // copy result from DEVICE to HOST.
  cudaMemcpy(&count, dev_count, sizeof(int), cudaMemcpyDeviceToHost);
  cudaMemcpy(indices, dev_indices, count*sizeof(uint2), cudaMemcpyDeviceToHost);
//*/
  for (int i = 0; i < count; ++i) {
    cout << left 
         << setw(3) << indices[i].x << ":"
         << setw(10) << brecords[indices[i].y].height << " , "
         << setw(3) << indices[i].y << ":"
         << setw(10) << arecords[indices[i].x].weight << endl;
  }
//*/
  cout << count << " records found." << endl;

  cudaFree(dev_arecords);
  cudaFree(dev_brecords);
  cudaFree(dev_count);
  cudaFree(dev_indices);
  cudaDeviceReset();

}

f:id:Episteme:20170201191910p:plain

Device-memoryの制限がよりいっそうキツくなります。Device-memoryに載せる配列がふたつになりますし、検索結果(intのペア)の格納域がNa*Nb に比例して増えますもん。

on-memory database の検索 : SELECT WHERE ~

週末の退屈しのぎにコード書き。

struct record {
  float height;
  float weight;
};

なんてな struct の配列 records[N] から weight < 0.1 を満たすものを列挙します。
SQLでいうところの

SELECT * 
  FROM records
  WHERE weight < 0.1;

ですわな。
これをCUDAでやるとなると records[N] および 検索結果を格納するindex列 int indices[N] はDevice上に置くことになります。
WHERE 節に与える条件式は device-lambda 使って楽しました。

// nvcc fake_DB.cpp --expt-extended-lambda 
// CUDA 8.0

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

struct record {
  float height;
  float weight;
};

template<typename Where>
__global__ void kernel_select(
    const record* records, unsigned int size,
    Where condition, int* count, int* indices) {

  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if (i == 0) *count = 0; // 一回だけ、 *count を 0-clear
  __syncthreads();

  while (i < size) {
    // i番目のrecordが WHERE条件を満たしているなら
    // indicesにiを追加
    if ( condition(records[i]) ) {
      indices[atomicAdd(count, 1)] = i;
    }
    i += gridDim.x * blockDim.x;
  }
}

#include <random>
#include <iostream>
#include <iomanip>
#include <algorithm>
using namespace std;

int main() {
  const int N = 100;
  record records[N];

  // fill records[] with random numbers
  mt19937 gen;
  uniform_real_distribution<float> dist;
  generate_n(records, N, [&]() { return record{ dist(gen),dist(gen) }; });

  // result
  int indices[N];
  int count;

  // allocate device-mem.
  record* dev_records;
  cudaMalloc(&dev_records, N * sizeof(record));
  int* dev_indices;
  cudaMalloc(&dev_indices, N * sizeof(int));
  int* dev_count;
  cudaMalloc(&dev_count, sizeof(int));

  // copy records from HOST to DEVICE
  cudaMemcpy(dev_records, records, N * sizeof(record), cudaMemcpyHostToDevice);

  // SELECT * FROM dev_records
  // WHERE weight < 0.1
  auto where = [] __device__(const record& rec) { return rec.weight < 0.1f; };

  dim3 grid = (N + 255) / 256;
  dim3 block = 256;
  kernel_select <<<grid, block >> >(
      dev_records, N,
      where,
      dev_count, dev_indices);

  // copy result from DEVICE to HOST.
  cudaMemcpy(&count, dev_count, sizeof(int), cudaMemcpyDeviceToHost);
  cudaMemcpy(indices, dev_indices, count * sizeof(int), cudaMemcpyDeviceToHost);
//*/
  for (int i = 0; i < count; ++i) {
    cout << left 
         << setw(10) << records[indices[i]].height << " : "
         << setw(10) << records[indices[i]].weight << endl;
  }
//*/
  cout << count << " records found." << endl;

  cudaFree(dev_records);
  cudaFree(dev_count);
  cudaFree(dev_indices);
  cudaDeviceReset();

}

f:id:Episteme:20170130214453p:plain

メモリの乏しいDeviceでは扱えるrecord数をそんなに大きくはできません。まぢに使うとなればrecords[N]を分割し、Deviceに乗せては検索を繰り返すことになるでしょうね。

Grid-stride loop

NVIDIAの開発系BlogになるほどなーなTips: Grid-stride loop を見つけたので軽く解説。 (元ネタはコチラ)

毎度おなじみ SAXPY : ベクトルの積和演算 Y[i] = a * X[i] + Y[i] (i = 0,1 ... n-1) これをCPUでフツーに行うとき、iについてくるくる回すloopで実装しますわな。

void saxpy(int n, float a, const float* x, float* y) {
  for ( int i = 0; i < n; ++i ) {
    y[i] = a * x[i] + y[i];
  }
}

かたやこいつをGPU(CUDA)でやるときの常套手段は "iによるn回のloop" を "n個のスレッド" に置き換えることで コードからloopを取っ払います。

__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < n ) {
    y[i] = a * x[i] + y[i];
  }
}

void saxpy(int n, float a, const float* x, float* y) {
  kernel_saxpy<<<(n+255)/256, 256>>>(n, a, x, y);
}

Grid-stride loopは両者の中間みたいなやりかたで、for-loopを温存します。

__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
  for ( int i = blockDim.x * blockIdx.x + threadIdx.x;
            i < n;
            i += gridDim.x * blockDim.x ) {
    y[i] = a * x[i] + y[i];
  }
}

このとき iの初期値はスレッドの通し番号、loopごとの増分はスレッドの総数すなわち gridDim.x * blockDim.x です。するってーとスレッド総数が Tのとき、

 スレッド#0 : i = 0, T  , 2T ...  
 スレッド#1 : i = 1, T+1, 2T+1 ...  
 スレッド#2 : i = 2, T+2, 2T+2 ...  
 (以下同文)

に対して y[i] を計算してくれる、と。

こうしておくとなにがオイシイかというと、kernel-call時のスレッド数を自由に指定できます。grid * blockがnより小さくてもちゃんと長さnのベクトルについて計算してくれる。

なのでたとえばGPUが持ってるプロセッサ数に応じてgrid/block数をチューニングできます。

void saxpy(int n, float a, const float* x, float* y) {
  int numSMs;
  cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
  // SMあたり32*256=8192個のスレッドを起こして計算する
  kernel_saxpy<<<numSMs*32, 256>>>(n, a, x, y);
}

極端な話 kernel_saxpy<<<1, 1>>>(n, a, x, y); でも問題ありません。このときスレッド一個で処理しますから多数のスレッドが絡む問題を起こさないので本来得られるべき結果が出てくるわけで、デバッグ時の検証に重宝しますね。 加えて kernel内で'prinf'したとき出力がお行儀よく並びますからデバッグが楽ですよね。

おためしに書いてたコードがこちらになります:

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

__global__ void kernel_saxpy(int n, float a, const float* x, float* y) {
  for ( int i = blockDim.x * blockIdx.x + threadIdx.x;
            i < n;
            i += gridDim.x * blockDim.x ) {
    y[i] = a * x[i] + y[i];
  }
}

#include <iostream>
#include <random>
#include <algorithm>
#include <ctime>

int main() {
  using namespace std;

  const int N = 100000;
  
  float* x;
  float* y0;
  float* y;
  float  a;

  cudaMallocManaged(&x, N*sizeof(float));
  cudaMallocManaged(&y, N*sizeof(float));
  cudaMallocManaged(&y0, N*sizeof(float));

  mt19937 gen(time(nullptr));
  uniform_real_distribution<float> dist;
  auto rand = [&]() { return dist(gen); };
  generate_n(x, N, rand);
  generate_n(y, N, rand);
  copy(y, y+N, y0);
  a = rand();

  // 単一スレッドで計算する
  kernel_saxpy<<<1,1>>>(N, a, x, y0);

  int numSMs;
  cudaDeviceGetAttribute(&numSMs, cudaDevAttrMultiProcessorCount, 0);
  cout << numSMs << " SMs" << endl;
  // SMあたり32*256=8192個のスレッドを起こして計算する
  kernel_saxpy<<<numSMs*32, 256>>>(N, a, x, y);
  cudaDeviceSynchronize();

  if ( equal(y0, y0+N, y) ) {
    cout << "ok." << endl;
  }

  cudaFree(x);
  cudaFree(y);
  cudaFree(y0);
  cudaDeviceReset();

}

転置しながらコピーする

cuBLAS が扱う行列要するに二次元配列のメモリレイアウトはFortranC/C++では大きく異なります。 行列

│ 1 2 3 │
│ 4 5 6 │

がメモリ上にリニアに展開されたとき、
Fortranでは 1 4 2 5 3 6 の順、上から下/左から右に並ぶのに対し、
C/C++ では 1 2 3 4 5 6 の順、左から右/上から下に並びます。

上から下(行:row方向)が左から右(列:column方向)より優先するFortran流レイアウトをrow-major、反対のC/C++流レイアウトをcolumn-majorっていいます。

f:id:Episteme:20161218131439j:plain

cuBLASの基となった本家BLASはもともとFortranで実装されたライブラリであるがため、行列表現はFortranに倣ってrow-majorになってます。したがって上記の行列をC/C++の二次元配列で表現すると、

  //      列 行
  float A[3][2] = {
    { 1.0f, 4.0f },
    { 2.0f, 5.0f },
    { 3.0f, 6.0f }
  };

となります。行列を転置(行/列を交換)せにゃならんのです。

これがどうにも居心地よろしくありません。Host上ではcolumn-major、cuBLASが扱うDevice上ではrow-majorで表現し、Host-Device間のコピーついでに転置する transpose_copy をこしらえてみました。

#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <iostream>

template<typename T>
cudaError_t transpose_copy(       T* dst, int ldDst, // コピー元
                            const T* src, int ldSrc, // コピー元
                                  int m,  int n,     // 行, 列
                            cudaStream_t stream = nullptr) {
  cudaError_t status = cudaSuccess;
  for ( int i = 0; i < n && status == cudaSuccess; ++i ) {
    status = cudaMemcpy2DAsync(dst, sizeof(T), 
                               src, sizeof(T)*ldSrc, 
                               sizeof(T), m, 
                               cudaMemcpyDefault, stream);
    ++src;
    dst += ldDst;
  }
  return status;
}

int main() {
  const int LD = 8;

  // 2行3列のA と
  float A[LD][LD] = {
    { 1.0f, 1.1f, 1.2f },
    { 1.3f, 1.4f, 1.5f }
  };
  // 3行2列のB の積を
  float B[LD][LD] = {
    { 2.0f, 2.1f },
    { 2.2f, 2.3f },
    { 2.4f, 2.5f }
  };
  // 2行2列のC に求める
  float C[LD][LD];

  float* dA;
  float* dB;
  float* dC;
  cudaMalloc(&dA, LD*LD*sizeof(float));
  cudaMalloc(&dB, LD*LD*sizeof(float));
  cudaMalloc(&dC, LD*LD*sizeof(float));

  // Host->Device
  transpose_copy(dA, LD, (float*)A, LD, 2, 3);
  transpose_copy(dB, LD, (float*)B, LD, 3, 2);
  // Cは0で埋めておく
  cudaMemset(dC, 0, LD*LD*sizeof(float));

  cublasHandle_t handle;
  cublasCreate(&handle);

  // GEMM: C = αAB + βC (α=1, β=0 なら C = AB)
  float alpha = 1.0f;
  float beta  = 0.0f;
  cublasSgemm(handle,
              CUBLAS_OP_N, CUBLAS_OP_N,
              2, 2, 3,
              &alpha, dA, LD, dB, LD, 
              &beta, dC, LD );
  // Device->Host
  transpose_copy((float*)C, LD, dC, LD, 2, 2);
  cudaDeviceSynchronize();

  for ( int y = 0; y < 2; ++y ) {
    for ( int x = 0; x < 2; ++x ) {
      std::cout << C[y][x] << ' ';
    }
    std::cout << std::endl;
  }

  cublasDestroy(handle);
  cudaFree(dA);
  cudaFree(dB);
  cudaFree(dC);

}

f:id:Episteme:20161218130218p:plain

nvgraphPagerankでPagerankを求める

CodeZinePagerankアルゴリズムの解説記事を書きました。
推移確率行列:H に対し

 λI = HI

を満たす固有ベクトル(eigenvector):I がすなわち Pagerank となります。
※ Iは確率ベクトルなので総和は常に1、なので固有値(eigenvalue):λ は1.0 です。

この問題を解いてくれるのが nvGRAPH の nvgraphPagerank、CodeZineで使ったサンプルデータで nvgraphPagerank を使ってみましたよ。
※ codeの半分は COO-format から CSC-format への変換でできてます。

f:id:Episteme:20161206203253p:plain

#include <cuda_runtime.h>
#include <nvgraph.h>

#include <iostream>
#include <algorithm>

using namespace std;

void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float*bookmark, float* probs, int iter);

template<typename T>
void dump(T* data, int n) {
  T* items = new T[n];
  cudaMemcpy(items, data, n*sizeof(T), cudaMemcpyDefault);
  for_each( items, items+n, [](T item) { cerr << item << " "; });
  cerr << endl;
  delete[] items;
}

int main() {
  const float r0 = 0.0f;
  const float r1 = 1.0f;
  const float r2 = 1.0f / 2.0f;
  const float r3 = 1.0f / 3.0f;
  const float r8 = 1.0f / 8.0f;

  // 元ネタ: COO-format
  const int nvertices =  8; // ページ数
  const int nedges    = 17; // リンク数
  int source_indices[nedges] = // リンク元
   {  0,  0,  1,  2,  2,  3,  3,  3,  4,  4,  4,  5,  6,  6,  6,  7,  7 }; 
  int destination_indices[nedges] = // リンク先
   {  1,  2,  3,  1,  4,  1,  4,  5,  5,  6,  7,  7,  0,  4,  7,  5,  6 };
  float weights[nedges] = // 各リンクの重み
   { r2, r2, r1, r2, r2, r3, r3, r3, r3, r3, r3, r1, r3, r3, r3, r2, r2 };

  // リンク先が少なくともひとつあるページは 1.0f、
  // リンク先のない行き止まりページは 0.0f
  float bookmark[nedges];
  fill_n(bookmark, nedges, 1.0f);
  for ( int dest : destination_indices ) {
    bookmark[dest] = 0.0f;
  }

  // 上記の元ネタを coo_topology に転写
  nvgraphCOOTopology32I_st coo_topology;
  float* coo_weights;
  cudaMalloc(&coo_topology.source_indices     , nedges*sizeof(int));
  cudaMalloc(&coo_topology.destination_indices, nedges*sizeof(int));
  cudaMalloc(&coo_weights, nedges*sizeof(float));

  coo_topology.nvertices = nvertices;
  coo_topology.nedges    = nedges;
  cudaMemcpy(coo_topology.source_indices     , source_indices     , nedges*sizeof(int), cudaMemcpyDefault);
  cudaMemcpy(coo_topology.destination_indices, destination_indices, nedges*sizeof(int), cudaMemcpyDefault);
  coo_topology.tag = NVGRAPH_UNSORTED;
  cudaMemcpy(coo_weights, weights, nedges*sizeof(float), cudaMemcpyDefault);

  // nvGRAPH の出番はここから
  nvgraphHandle_t handle;
  nvgraphCreate(&handle);
  cudaDataType_t data_type = CUDA_R_32F; // 32bit-real (=float)

  {
  // Pagerankを求めるには CSR-formatに変換せにゃならん
  nvgraphCSCTopology32I_st csc_topology;
  float* csc_weights;
  cudaMalloc(&csc_topology.source_indices, nedges*sizeof(int));
  cudaMalloc(&csc_topology.destination_offsets, (nvertices+1)*sizeof(int));
  cudaMalloc(&csc_weights, nedges*sizeof(float));

  nvgraphConvertTopology(handle, 
                         NVGRAPH_COO_32, &coo_topology, coo_weights,  // in  type, topology, weigts
                         &data_type,                                  // type of weights
                         NVGRAPH_CSC_32, &csc_topology, csc_weights); // out type, topology, weigts

  cerr << "CSC : " <<endl;
  cerr << "nvertices           = " << csc_topology.nvertices << endl;
  cerr << "nedges              = " << csc_topology.nedges    << endl;
  cerr << "source_indices      = "; dump(csc_topology.source_indices     , csc_topology.nedges);
  cerr << "destination_offsets = "; dump(csc_topology.destination_offsets, csc_topology.nvertices+1);
  cerr << "weights             = "; dump(csc_weights, csc_topology.nedges);

  // Pagerankを求めて出力
  float probs[nvertices] = { r1, r0, r0, r0, r0, r0, r0, r0 };
  pagerank(handle, &csc_topology, csc_weights, bookmark, probs, 100);
  cout << "\nPagerank: \n";
  for ( int i = 0; i < nvertices; ++i ) {
    cout << "page[" << i << "] : " << probs[i] << endl;
  }

  cudaFree(csc_topology.source_indices);
  cudaFree(csc_topology.destination_offsets);
  cudaFree(csc_weights);
  }

  cudaFree(coo_topology.source_indices);
  cudaFree(coo_topology.destination_indices);
  cudaFree(coo_weights);

  nvgraphDestroy(handle);
  cudaDeviceReset();
}

void pagerank(nvgraphHandle_t handle, nvgraphCSCTopology32I_t topology, float* weights, float* bookmark, float* probs, int iter) {
  nvgraphGraphDescr_t descr;
  // setup
  nvgraphCreateGraphDescr(handle, &descr);

  // topology
  nvgraphSetGraphStructure(handle, descr, topology, NVGRAPH_CSC_32);
  cudaDataType_t dimT[2] = { CUDA_R_32F, CUDA_R_32F };
  nvgraphAllocateVertexData(handle, descr, 2, dimT);
  nvgraphAllocateEdgeData(handle, descr, 1, dimT);
  // weight
  nvgraphSetEdgeData(handle, descr, weights, 0);
  // bookmark
  nvgraphSetVertexData(handle, descr, bookmark, 0);
  // initial guess
  nvgraphSetVertexData(handle, descr, probs, 1);

  // pagerank!!
  float damping_factor = 0.85f; // 85%はリンクを辿り,残る15%は気ままにジャンプ
  nvgraphPagerank(handle, descr,
    0,       // weight    : EdgeData[0] 
    &damping_factor, 
    0,       // bookmark  : VertexData[0]
    1,       // has guess
    1,       // pagerank  : VertexData[1]
    1.0e-6f, // tolerance 誤差 : これより小さくなったら計算終了
    iter     // 最大繰り返し回数
  );
  // get result
  nvgraphGetVertexData(handle, descr, probs, 1);

  // teardown
  nvgraphDestroyGraphDescr(handle, descr);
}

f:id:Episteme:20161206203340p:plain

thrust::for_eachによるヒストグラム

Thrust、個人的にはかなり気に入ってるライブラリなんだけど、なんだか影が薄いんですよね。 cuFFTやcuDNNなどなど、僕らにはとても書けそうにないガチにヘビィな連中とは異なり、Thrustが提供するのは copy/search/sortといった基本的なアルゴリズム群だからでしょうか。

CUDAが、というかnvccが device lambda をサポートしてくれたおかげで、Thrustの使い心地がぐっと良くなりました。その一例。

ヒストグラム(度数分布表)を作ります。たくさんの数値データがあって、それらのとる値の範囲を等間隔に刻んだスロット列を用意します。各数値をスロットに振り分けたとき、各スロットにいくつずつ納まるかってやつ。

アルゴリズムはこんなカンジ:

  float score[N]; // 数値データ列(input)
  int   hist[S]; // ヒストグラム(output) 初期値はall-0
  for ( score 内の各数値 val について ) {
    if ( valがスロットiの範囲にあるなら } {
      ++hist[i};
    }
  }

これをkernelで実装するのはちっとも難しくありません。けどもたったこれだけのためにわざわざkernel書くのが面倒なわけで、Thrustでさらっと実装してみました。

#include <thrust/device_vector.h>
#include <thrust/for_each.h>

#include <curand.h>

#include <iterator>
#include <iostream>
#include <ctime>

int main() {
  using namespace std;

  const int N = 100000;
  thrust::device_vector<float> dscore(N);
  thrust::device_vector<int>   dhist(101, 0);

  { // 平均 50.0, 標準偏差 15.0 の正規分布乱数で dscoreを埋める 
    curandGenerator_t gen;
    curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT);
    curandSetPseudoRandomGeneratorSeed(gen, static_cast<unsigned long long>(time(nullptr)));
    curandGenerateNormal(gen, dscore.data().get(), N, 50.30f, 15.0f);
    curandDestroyGenerator(gen);
  }

  // dscoreの各値を四捨五入し、0以上/100以下ならばヒストグラムに加える
  int* histPtr = dhist.data().get(); // for_each内のlambdaがコレをキャプチャする
  thrust::for_each(begin(dscore), end(dscore),
                   [=] __device__ (float val) -> void {
                     int i = static_cast<int>(val+0.5f); // 四捨五入して
                     if ( i >= 0 && i <= 100 ) { // 0以上100以下なら、当該スロットを+1する
                       atomicAdd(histPtr+i, 1); // ++histPtr[i] だと data-race起こすよきっと!
                     }
                   });
//*/
  for ( int item : dhist ) {
    cout << item << endl;
  }
//*/
}

thrust::for_each(こっから, ここまで, なんかする) の"なんかする"を devie lambda で与えています。ヒストグラムを求める領域:dhistのポインタをlambda式の値キャプチャで食わせています。

※ device lambdaでキャプチャできるのは値キャプチャだけ。参照キャプチャするとコンパイル・エラーとなります。Hostにある変数の参照(つまりポインタ)をDevice側に持ち込んでも意味ないもんね。

結果をCSVにリダイレクトし、excelに食わせたものがこちらになります。

f:id:Episteme:20161130192402p:plain

zip_iteratorによるSoAのソート

ふたつの配列: int x[N], y[N] があって、このふたつを連動してソートしたい、ソート順は「xの小さい順だけど、xが同じときはyの小さい順」...あるある。 CUDA世界ではデータを SoA(Structure of Array)で構成することが多いのでなおさらよくあるシチュエーションです。

thrust::sort()はひとつの配列をソートするのは簡単だけど、複数の配列を連動してソートさせるにはどーすりゃいいんだと。

thrust::zip_iteratorっていう小賢しいiteratorが用意されています。 複数のiteratorを一本のzip_iteratorにまとめることができ、zip_iteratorの指す先が移動するとまとめられたiteratorそれぞれが同じだけ移動します。 加えてzip_iteratorに対しoperator*()で値を取り出すと各iteratorの指す先の値がtupleにまとめられて返ってくるですよ。

#include <thrust/iterator/zip_iterator.h>
#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/sort.h>

#include <random>
#include <algorithm>
#include <numeric>

#include <iterator>
#include <iostream>

int main() {
  using namespace std;

  const int N = 90;
  thrust::host_vector<int> x(N);
  thrust::host_vector<int> y(N);

  { // テキトーな組を生成する
    mt19937 gen;
    iota(begin(x), begin(x)+N/2, 10);
    iota(begin(x)+N/2, end(x), 10);
    iota(begin(y), end(y), 10);
    shuffle(begin(x), end(x), gen);
    shuffle(begin(y), end(y), gen);
  }

  cout << "--- before:\n";
  for ( int i = 0; i < N; ++i ) { cout << x[i] << '-' << y[i] << ' '; }
  cout << endl;

  // deviceにコピー
  thrust::device_vector<int> dx = x;
  thrust::device_vector<int> dy = y;

  // dx, dy のbegin/end() をzip_iteratorでまとめ、
  auto first = thrust::make_zip_iterator(thrust::make_tuple(begin(dx), begin(dy)));
  auto last  = thrust::make_zip_iterator(thrust::make_tuple(end(dx)  , end(dy)  ));

  // dx,dy を連動してソートする
  thrust::sort(first, last); 

  // hostに書き戻し
  x = dx;
  y = dy;
  cout << "--- after(ascending):\n";
  for ( int i = 0; i < N; ++i ) { cout << x[i] << '-' << y[i] << ' '; }
  cout << endl;

  // もう一度、今度は(比較ファンクタを与えて)降順で。
  thrust::sort(first, last, 
               [] __device__ (const auto& a, const auto& b) -> bool { return b < a; });

  x = dx;
  y = dy;
  cout << "--- after(descending):\n";
  for ( int i = 0; i < N; ++i ) { cout << x[i] << '-' << y[i] << ' '; }
  cout << endl;
}

f:id:Episteme:20161128195150p:plain