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

東方算程譚

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

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();

}

転置しながらコピーする

Memory

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を求める

nvGRAPH

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

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のソート

Thrust

ふたつの配列: 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

Thrust で 生ポを扱うには

Thrust

Windows/Visual Studio版ではおなじみの「配列の足し算」: c[i] = a[i] + b[i] (i = 0, 1, ...) を Thrust使って書いてみました。

/*
  DO NOT FORGET: --expt-extended-lambda option to nvcc
 */

#include <cstdio>
#include <cassert>

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <thrust/transform.h>

void addWithThrust(thrust::host_vector<int>& c, 
                   const thrust::host_vector<int>& a, 
                   const thrust::host_vector<int>& b);

int main() {
  using namespace std;

  const int arraySize = 5;
  const int a_data[arraySize] = {  1,  2,  3,  4,  5 };
  const int b_data[arraySize] = { 10, 20, 30, 40, 50 };

  thrust::host_vector<int> a(a_data, a_data + arraySize);
  thrust::host_vector<int> b(b_data, b_data + arraySize);
  thrust::host_vector<int> c(arraySize);

  // Add vectors in parallel.
  addWithThrust(c, a, b);

  printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
        c[0], c[1], c[2], c[3], c[4]);

  cudaDeviceReset();
}

void addWithThrust(      thrust::host_vector<int>& c, 
                   const thrust::host_vector<int>& a, 
                   const thrust::host_vector<int>& b) {
  using namespace std;

  assert( a.size() == b.size() );
  assert( a.size() == c.size() );

  thrust::device_vector<int> dev_a = a; // copy Host to Device
  thrust::device_vector<int> dev_b = b; // copy Host to Device
  thrust::device_vector<int> dev_c(c.size());

  thrust::transform(begin(dev_a), end(dev_a), // dev_a for input
                    begin(dev_b),             // dev_b for input
                    begin(dev_c),             // dev_c for output
                    [] __device__ (int x, int y) -> int { return x + y; }); // z = x + y

  c = dev_c; // copy Device to Host
}

楽ちんですねー、thrust::host_vector/device_vector にメモリ管理を任せてしまえば cudaMallocはコンストラクタ、cudaMemcpyoperator=()がやってくれるのでコードが実に涼しげです。

Thrustのコードを追いかけてみたところ、thrust:host_vectorのメモリ確保/解放は new/deletethrust::device_vectorcudaMalloc/cudaFreeが使われてるようです。

ここまで簡単になるんだからHost/Deviceで共用できるManaged Memory使えばもっと楽できんじゃないかと試してみたですよ。

#include <cstdio>
#include <cassert>

#include <algorithm>

#include <cuda_runtime.h>

#include <thrust/transform.h>

void addWithThrust(int* c,
                   const int* a,
                   const int* b,
                   int size);

int main() {
  using namespace std;

  const int arraySize = 5;
  const int a_data[arraySize] = {  1,  2,  3,  4,  5 };
  const int b_data[arraySize] = { 10, 20, 30, 40, 50 };

  int* a; cudaMallocManaged(&a, arraySize*sizeof(int)); std::copy( a_data, a_data + arraySize, a);
  int* b; cudaMallocManaged(&b, arraySize*sizeof(int)); std::copy( b_data, b_data + arraySize, b);
  int* c; cudaMallocManaged(&c, arraySize*sizeof(int)); std::fill( c, c+arraySize, 0);

  // Add vectors in parallel.
  addWithThrust(c, a, b, arraySize);

  printf("{1,2,3,4,5} + {10,20,30,40,50} = {%d,%d,%d,%d,%d}\n",
        c[0], c[1], c[2], c[3], c[4]);

  cudaFree(a);
  cudaFree(b);
  cudaFree(c);

  cudaDeviceReset();
}

void addWithThrust(int* c,
                   const int* a,
                   const int* b,
                   int size) {

  thrust::transform(a, a+size, // a for input
                    b,         // b for input
                    c,         // c for output
                    [] __device__ (int x, int y) -> int { return x + y; }); // z = x + y
  cudaDeviceSynchronize();

}

f:id:Episteme:20161121192437p:plain

あらら、ぜーんぜん動いてません。原因はココ:

  thrust::transform(a, a+size, // a for input
                    b,         // b for input
                    c,         // c for output
                    [] __device__ (int x, int y) -> int { return x + y; }); // z = x + y

Thrustが提供する多くの関数は、引き渡されるイテレータの型に応じてHostで行うかDeviceで行うかを静的に判断しています。このコードではイテレータとして生のポインタを渡してるんですけど、Thrustではイテレータが(生の)ポインタであったならHostで実行すべきものと判断してるっポいです。

解決策はふたつ。ひとつは第一引数に実行ポリシーを明示的に与えること。

#include <thrust/execution_policy.h>
...
void addWithThrust(int* c,
                   const int* a,
                   const int* b,
                   int size) {

  thrust::transform(thrust::device,
                    a, a+size, // a for input
                    b,         // b for input
                    c,         // c for output
                    [] __host__ __device__ (int x, int y) -> int { return x + y; }); // z = x + y
  cudaDeviceSynchronize();
}

もうひとつは、生のポインタを thrust::device_ptr に変換すること。

#include <thrust/device_ptr.h>
...
void addWithThrust(int* c,
                   const int* a,
                   const int* b,
                   int size) {

  thrust::device_ptr<int> pa(const_cast<int*>(a));
  thrust::device_ptr<int> pb(const_cast<int*>(b));
  thrust::device_ptr<int> pc(c);
  thrust::transform(pa, pa+size, // a for input
                    pb,          // b for input
                    pc,          // c for output
                    [] __host__ __device__ (int x, int y) -> int { return x + y; }); // z = x + y
  cudaDeviceSynchronize();

}

COO-format から CSR-format への変換

cuSPARSE

m行n列の行列、成分数は m*n個。m,nが大きくなるとメモリの使用量がシャレになりません。けどもほとんどの成分が0であるなら、0成分を省略することでメモリ消費をぐっと抑えることができます。

ほとんどの成分が0のスカスカな行列を疎行列(Sparse matrix)と称します。その疎行列を対象にした各種演算をCUDAでやらせるライブラリがcuSPARSEです。

疎行列を表現するには (行index, 列index, 成分) の組を非0成分数だけ並べればよいですな。SoA(Structure of Array)なレイアウトだと

  int   rowInd[nnz]; // 行index
  int   colInd[nnz]; // 列index
  float val[nnz];    // 成分

こんなカンジ、nnzは非0成分の数(number of nonzero)。この形式をCOO-format(Coodinate format)っていいます。

くSPARSEが提供する演算には疎行列を引数に与えるんですけど、上記COO-formatのままでは引数に与えることができず、CSR-format(Compressed Sparse Row format)に変換せにゃなりません。その方法。

まずCOO-formatの各組を行(row),列(col)の順に昇順でソートします。行(row)の小さい順、行が同じなら列(col)の小さい順。当然成分(val)も連動して入れ替えを行います。

こんな元データが

f:id:Episteme:20161116191559p:plain

このようにソートされます。

f:id:Episteme:20161116191605p:plain

これが row-major でソートされた COO-format。

行indexは同じ値が連続しますね、row-majorにソートしたんだから当然ですけど。 同じ値が連続するってことはすなわち冗長であり圧縮できるってことです。

rowInd[nnz] を rowPtr[m+1] に圧縮した結果がコチラ。

f:id:Episteme:20161116191611p:plain

これがCSR-format(Compressed Sparse matrix Row format)。

COO-format から CSR-format への変換関数群は cuSPARSE の中に用意されています。

#include <cuda_runtime.h>
#include <cusparse.h>
#include <algorithm>
#include <numeric>
#include <random>
#include <iostream>
#include <iomanip>

using namespace std;

int main() {
  const int m   = 10; // 行数
  const int n   = 20; // 列数
  const int nnz = 30; // 非0成分数(numner of non-zero)

  // 疎行列のCOO(Coodinate format)
  int*   rowInd;
  int*   colInd;
  float* cooVal;
  cudaMallocManaged(&rowInd, nnz*sizeof(int));
  cudaMallocManaged(&colInd, nnz*sizeof(int));
  cudaMallocManaged(&cooVal, nnz*sizeof(float));

  { // 乱数で疎行列を生成する
    int inx[m*n];
    iota(begin(inx), end(inx), 0);
    shuffle(begin(inx), end(inx), mt19937());
    for ( int i = 0; i < nnz; ++i ) {
      rowInd[i] = inx[i]/n;
      colInd[i] = inx[i]%n;
      cooVal[i] = rowInd[i] + colInd[i]*0.01f;
    }
  }

  cout << "COO-format\n"
          "# ( rowInd, colInd, val)\n";
  for ( int i = 0; i< nnz; ++i ) {
    cout << setw(3) << i << " ( "
         << setw(3) << rowInd[i] << ',' 
         << setw(3) << colInd[i] << ", " 
         << cooVal[i] 
         << " )" << endl;
  }
  cout << endl;

  // sort後の成分領域を確保
  float* csrVal;
  cudaMallocManaged(&csrVal, nnz*sizeof(float));

  // rowPtr(圧縮されたrowInd)領域を確保 大きさは m+1
  int* rowPtr;
  cudaMallocManaged(&rowPtr, (m+1)*sizeof(int));

  // cuSPARSE 初期化
  cusparseHandle_t handle;
  cusparseCreate(&handle);

  // coosortに必要なバッファの大きさを求め、
  size_t bufferSize;
  cusparseXcoosort_bufferSizeExt(handle, m, n, nnz, rowInd, colInd, &bufferSize);
  // バッファを確保する
  void* sortBuffer;
  cudaMalloc(&sortBuffer, bufferSize);

  // permutation領域を確保し、
  int* permutation;
  cudaMalloc(&permutation, nnz*sizeof(int));
  // 0, 1, 2,... で埋める
  cusparseCreateIdentityPermutation(handle, nnz, permutation);

  // 行、列の順で昇順にソートし、
  cusparseXcoosortByRow(handle, m, n, nnz, rowInd, colInd, permutation, sortBuffer);
  cudaFree(sortBuffer);

  // 成分がソート順に対応するようpermutationを基に配置する
  cusparseSgthr(handle, nnz, cooVal, csrVal, permutation, CUSPARSE_INDEX_BASE_ZERO);
  cudaFree(permutation);

  // rowIndを圧縮してrwoPtrに求める
  cusparseXcoo2csr(handle, rowInd, nnz, m, rowPtr, CUSPARSE_INDEX_BASE_ZERO);
  cudaDeviceSynchronize();

  cout << "COO-format row-major sorted\n"
          "# ( rowInd, colInd, val)\n";
  for ( int i = 0; i< nnz; ++i ) {
    cout << setw(3) << i << " ( "
         << setw(3) << rowInd[i] << ',' 
         << setw(3) << colInd[i] << ", " 
         << csrVal[i] 
         << " )" << endl;
  }
  cout << endl;

  cout << "CSR-format\n"
          "# ( rowPtr, colInd, val)\n";
  for ( int i = 0; i< nnz; ++i ) {
    cout << setw(3) << i; ;
    if ( i <= m )
      cout << " ( " << setw(3) << rowPtr[i] << ',';
    else
      cout << "     ( ";
    cout << setw(3) << colInd[i] << ", " 
         << csrVal[i] 
         << " )" << endl;
  }
  cout << endl;

  cudaFree(rowInd);
  cudaFree(rowPtr);
  cudaFree(colInd);
  cudaFree(cooVal);
  cudaFree(csrVal);
  cusparseDestroy(handle);
  cudaDeviceReset();

}