ニックネーム:επιστημη
都道府県:Kanagawa, JAPAN

»くわしく見る
 
multithread環境下の default-stream (そのに)
2016年07月25日(月)
前回のつづきというか補足。

コンパイル・オプション --default-stream=per-thread を与えておくと、ホストthreadごとにdefault-streamが作られます。

threadごとに個別のstreamがあてがわれるなら、ホスト⇔デバイス間のメモリ・コピーとkernel実行がoverlapできるはずよね。ってことで実験。

ホスト側メモリを cudaHostAlloc で確保し、cudaMemcpyをcudaMemcpyAsyncに差し替え、二つのホストthreadでメモリ・コピーとkernel実行を繰り返す様子を nvvp で観察しました。



茶色の帯がメモリ・コピー、赤/青の帯がkernel実行です。
両者がoverlapし、処理時間的にはほぼkernel実行にかかる時間の分だけ速くなってる勘定になりました。

複数のホストthreadがCUDAを使うときの(streamを明示的に用意しないお手軽な)高速化テクニックとして使えそうです。

2016-07-25 15:48 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/122/

multithread環境下の default-stream
2016年07月19日(火)
CUDAにおけるstreamは受付窓口に並ぶ待ち行列みたいなもんで、窓口はひとつしかないのですがそこに並ぶ列は複数個作れます。

cudaMemcpyやkernel-call時にstreamを指定しない場合には、あらかじめ用意されているdefault-streamに並びます。

そこで実験、複数のHost-threadからCUDAを使うとどーなるんでしょ。
おなじみ cudaMemcpy(HtoD)/kernel-call/cudaMemcpy(DtoH)の一連のシーケンスをふたつのスレッドで同時進行させ、実行の様子をNVVPで観察します。
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

template<typename T, typename UnaryFunction>
__global__ void kernel_transform(T* y, const T* x, unsigned int size, UnaryFunction f) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) { y[i] = f(x[i]); }
}

template<typename T>
void launch_sincos(T* y, const T* x, unsigned int size, bool is_sin) {
  T* dx; cudaMalloc(&dx, size*sizeof(T));
  T* dy; cudaMalloc(&dy, size*sizeof(T));

  cudaMemcpy(dx, x, size*sizeof(int), cudaMemcpyHostToDevice);
  if ( is_sin )
    kernel_transform<<<(size+255)/256,256>>>(dy, dx, size, 
                        [] __device__ (T arg) { return sin(arg);} );
  else
    kernel_transform<<<(size+255)/256,256>>>(dy, dx, size, 
                        [] __device__ (T arg) { return cos(arg);} );
  cudaMemcpy(y, dy, size*sizeof(int), cudaMemcpyDeviceToHost);

  cudaFree(dx);
  cudaFree(dy);
}

#include <thread>
#include <vector>
#include <numeric>
#include <functional>

int main() {
  using namespace std;
  const unsigned int size = 100000;

  vector<float> x0(size);
  vector<float> x1(size);

  iota(begin(x0), end(x0), 0.0f);
  iota(begin(x1), end(x1), 0.0f);

  vector<float> y0(size);
  vector<float> y1(size);

  thread thr0(launch_sincos<float>, y0.data(), x0.data(), size, true );
  thread thr1(launch_sincos<float>, y1.data(), x1.data(), size, false);
  thr0.join(); thr1.join();

}



ふむ、複数のスレッドで処理してもdefault-streamは一本です。

ここでnvccのオプションに --default-stream=per-thread を指定すると


スレッドごとにdefault-streamが用意されるんですわ。



2016-07-19 08:51 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/121/

NVTX: ユーザ定義のtimeline (そのに)
2016年07月11日(月)
前回のつづき。

nvtxRangeStart() と nvtxRangeEnd() ではさむことで NVVPのタイムチャート上に timeine を乗っけることができるんだけど、もうひとつ、nvtxRangePush() と nvtxRangePop() ではさむってのがあります。

こちらは timeline を"入れ子"にできるので、関数コールのネストなんかを表示するのに向いてます。再帰的なソートを実装し、timelneをひっぱってみました。
#include <iostream>
#include <random>
#include <algorithm>
#include <numeric>
#include <vector>
#include <string>
#include <nvToolsExt.h>

// data[first] 〜 data[last-1] をソートする
template<typename T>
void my_sort(T* data, size_t first, size_t last) {
  nvtxRangePush((std::to_string(first) + " - " + std::to_string(last)).c_str());
  if ( last - first < 100 ) {
    // 要素数が少ないときはフツーにソート
    std::sort(data+first, data+last);
  } else {
    // 前半/後半をソートして
    size_t mid = (first + last) >> 1;
    my_sort(data, first, mid);
    my_sort(data, mid, last );
    // しかるのちマージする
    std::inplace_merge(data+first, data+mid, data+last);
  }
  nvtxRangePop();
}

int main() {
  using namespace std;
  const int N = 500;
  vector<float> data(N);

  iota(data.begin(), data.end(), 0.0f);
  shuffle(data.begin(), data.end(), mt19937());
  my_sort(data.data(), 0, N);
  cout << (is_sorted(data.begin(), data.end()) ? "ok\n" : "ng\n");

}

2016-07-11 08:08 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / NVTX / nvvp |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/120/

NVTX: ユーザ定義のtimeline
2016年07月04日(月)
NVVP(NVIDIA Visual Profiler)はその名のとおり、実行モジュールのプロファイリングを行い、メモリ・コピーやkernel実行などの処理時間をタイムチャートで見せてくれる便利ツールです。

んで、NVVPにはNVTX(NVIDIA Tools Extension)っていうおマケがついてます。コレ使うとユーザ定義のtimelineをタイムチャート上に乗っけてくれます。うれしいことにCUDAを使ってなくてもOK。

NVTXの在処は、
 NVTX_ROOT = C:\Program Files\NVIDIA Corporation\NvToolsExt
を基点として:
 include : <NVTX_ROOT>\include
 library : <NVTX_ROOT>\lib\x64\nvToolsExt64_1.lib
 bin(DLL): <NVTX_ROOT>\bin\x64\nvToolsExt64_1.dll
に配置されています。

vector<float> にデタラメな値をセットし / ソートして / 結果を検証する一連の処理をタイムチャートで図示してみます。

使い方は簡単、timelineの表示範囲を nvtxRangeStart() と nvtxRangeEnd() ではさむだけ。
#include <iostream>
#include <random>
#include <algorithm>
#include <numeric>
#include <vector>

#include <nvToolsExt.h>

int main() {
  using namespace std;
  const int N = 1000;
  vector<float> data(N);
  nvtxRangeId_t rid;

  rid = nvtxRangeStart("initialize"); // こっから
  iota(data.begin(), data.end(), 0.0f);
  shuffle(data.begin(), data.end(), mt19937());
  nvtxRangeEnd(rid); // ここまで

  rid = nvtxRangeStart("sort"); // こっから
  sort(data.begin(), data.end());
  nvtxRangeEnd(rid); // ここまで

  rid = nvtxRangeStart("certify"); // こっから
  cout << (is_sorted(data.begin(), data.end()) ? "ok\n" : "ng\n");
  nvtxRangeEnd(rid); // ここまで

}

コンパイルしてNVVP上で実行すると:


timelineをクリックすると所要時間も教えてくれます。


2016-07-04 09:05 | 記事へ | コメント(0) | トラックバック(0) |
| NVTX |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/119/

Thrust: key-value pair をソートする(そのさん)
2016年06月27日(月)
前々回前回のつづき。

今度はソート順を決定するkeyがひとつじゃないとき。
  int   key1[N];
  int   key2[N];
  float value[N];
こんな key1-key2-value を連動してソートします。key1昇順で並べ、key1が等しいならkey2昇順で並べます。
学年/組や苗字/名前でソートするってことね。

thrust::sort/sort_by_key はひとつのkeyで並べることしかできんけど、こいつらには比較関数オブジェクトを与えることができます。これを使ってですね...

まず int map[N] を用意し、0,1,2 ... N-1 で埋めておきます。んでもって thrust::sort でmap[N]をソートするんですけど、ここで与えた比較関数オブジェクトにはmap[N]の要素 a,b が与えられます。「a を b より手前に置きたいならtrueを返す」のが比較関数オブジェクトのルール。なので:

- key1[a] < key1[b] ならtrue
- key1[a] > key1[b] ならfalse
- key1[a] = key1[b] なら key2[a] < key2[b]

を返せば、key1をmajor-key, key2をminor-keyとしてソートしてくれます。
ソートが完了すればmap[N]は(前回と同様)要素の入れ替え表になってますから、これに従ってkey1, key2, value を入れ替えます。
/*
  DO NOT FORGET --expt-extended-lambda 
 */
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include <thrust/device_ptr.h>
#include <thrust/sort.h>

#include <iostream>
#include <vector>    // vector
#include <random>    // random
#include <utility>   // pair
#include <algorithm> // transform, shuffle

//#define GATHER_BY_THRUST

#ifdef GATHER_BY_THRUST
#include <thrust/gather.h>
#else
// out[i] = in[map[i]] for i = 0, 1, .. size-1
template<typename T>
__global__ void kernel_gather(const int* map, const T* in, T* out, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  int j;
  if ( i < size && (j = map[i]) < size ) {
    out[i] = in[j];
  }
}
#endif

// out[size] を 0, 1, 2, 3, ... で埋める
__global__ void kernel_iota(int* out, unsigned int size, int bias =0) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = i + bias;
  }
} 

int main() {
  using namespace std;

  const int n =2000;
  using record = pair<pair<int,int>,float>; // 2keys and value
  vector<record> aos(n);

  mt19937 gen;
  uniform_real_distribution<float> dist;
  // 元データをつくって
  for ( int i = 0; i < n; ++i ) { aos[i] = make_pair(make_pair(i/4,i%4), dist(gen)); }
  cout << "original:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << aos[i].first.first << '-' << aos[i].first.second << ':' 
         << aos[i].second << ' ';
  } 
  cout << "...\n";

  // かき混ぜる
  shuffle(aos.begin(), aos.end(), gen);
  cout << "shuffled:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << aos[i].first.first << '-' << aos[i].first.second << ':' 
         << aos[i].second << ' ';
  } 
  cout << "...\n";

  vector<int>   h_key1(n);
  vector<int>   h_key2(n);
  vector<float> h_value(n);

  // AOS から SOA にコピー
  for ( int i = 0; i < n; ++i ) {
    h_key1[i]  = aos[i].first.first;
    h_key2[i]  = aos[i].first.second;
    h_value[i] = aos[i].second;
  }

  int*   d_key1;
  int*   d_key2;
  float* d_value;

  void*  d_work;
  int*   d_map;

  cudaMalloc(&d_key1,  n*sizeof(int));
  cudaMalloc(&d_key2,  n*sizeof(int));
  cudaMalloc(&d_value, n*sizeof(float));

  cudaMalloc(&d_work,  n*max(sizeof(int),sizeof(float)));
  cudaMalloc(&d_map,   n*sizeof(int));

  // keyをdeviceにコピーし、mapを初期化
  cudaMemcpy(d_key1, h_key1.data(), n*sizeof(int), cudaMemcpyHostToDevice);
  cudaMemcpy(d_key2, h_key2.data(), n*sizeof(int), cudaMemcpyHostToDevice);
  kernel_iota<<<(n+255)/256,256>>>(d_map, n); // 0, 1, 2, ... n-1 をセット

  // key/mapをソートする
  thrust::device_ptr<int>   p_map(d_map);

  thrust::sort(p_map, p_map+n, 
              [=] __device__ (int a, int b) {
                return ( d_key1[a] < d_key1[b] ) ? true  : 
                      (( d_key1[b] < d_key1[a] ) ? false : 
                         d_key2[a] < d_key2[b]);
              });

  // keymapに基づいてkey1を入れ替え
  cudaMemcpy(d_work, d_key1, n*sizeof(float), cudaMemcpyDeviceToDevice);
#ifdef GATHER_BY_THRUST
  thrust::gather(p_map, p_map+n, 
                 thrust::device_ptr<int>((int*)d_work), 
                 thrust::device_ptr<int>(d_key1));
#else
  kernel_gather<<<(n + 255)/256, 256>>>(d_map, (int*)d_work, d_key1, n);
#endif

  // keymapに基づいてkey2を入れ替え
  cudaMemcpy(d_work, d_key2, n*sizeof(float), cudaMemcpyDeviceToDevice);
#ifdef GATHER_BY_THRUST
  thrust::gather(p_map, p_map+n, 
                 thrust::device_ptr<int>((int*)d_work), 
                 thrust::device_ptr<int>(d_key2));
#else
  kernel_gather<<<(n + 255)/256, 256>>>(d_map, (int*)d_work, d_key2, n);
#endif

  // keymapに基づいてvalueを入れ替え
  cudaMemcpy(d_work, h_value.data(), n*sizeof(float), cudaMemcpyHostToDevice);
#ifdef GATHER_BY_THRUST
  thrust::gather(p_map, p_map+n, 
                 thrust::device_ptr<float>((float*)d_work), 
                 thrust::device_ptr<float>(d_value));
#else
  kernel_gather<<<(n + 255)/256, 256>>>(d_map, (float*)d_work, d_value, n);
#endif

  // 結果をhostに書き戻して
  cudaMemcpy(h_key1.data(),  d_key1,  n*sizeof(int),   cudaMemcpyDeviceToHost);
  cudaMemcpy(h_key2.data(),  d_key2,  n*sizeof(int),   cudaMemcpyDeviceToHost);
  cudaMemcpy(h_value.data(), d_value, n*sizeof(float), cudaMemcpyDeviceToHost);

  // 確認
  cout << "sorted:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << h_key1[i] << '-' << h_key2[i] << ':' 
         << h_value[i] << ' ';
  } 
  cout << "...\n";

  // あとしまつ
  cudaFree(d_key1);
  cudaFree(d_key2);
  cudaFree(d_value);
  cudaFree(d_work);
  cudaFree(d_map);
  cudaDeviceReset();
}

2016-06-27 07:59 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/118/

Thrust: key-value pair をソートする(そのに)
2016年06月20日(月)
前回のつづき。thrust::sort_by_keyを使えばSOA(Structure of Array)なkey-value pairをソートすることができます。

では keyと連動して並び替えたいvalueがひとつじゃなかったら?
  int   key[N];
  float value1[N];
  float value2[N];
なんてな key-value1-value2 をソートするには?

...こんな方法はいかがでしょうか。

まず int map[N] を用意し、そのナカミを 0, 1, 2 ... N-1 で埋めておきます。しかるのち、thrust::sort_by_keyで key-map をソートします。

そーすっと、ソート完了後のmap[i]には"i番目にあるべき要素のindex"が入ってますわな。map[N]は要素の"入れ替え表"になってるわけですわ。

んでもって vaue1[N]とvalue2[N]の各要素をmap[N]の示すindexで移動すりゃいいじゃん、と。

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

#include <thrust/device_ptr.h>
#include <thrust/sort.h>

#include <iostream>
#include <vector>    // vector
#include <random>    // random
#include <utility>   // pair
#include <algorithm> // transform, shuffle

#ifdef GATHER_BY_THRUST
#include <thrust/gather.h>
#else
// out[i] = in[map[i]] for i = 0, 1, .. size-1
template<typename T>
__global__ void kernel_gather(const int* map, const T* in, T* out, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  int j;
  if ( i < size && (j = map[i]) < size ) {
    out[i] = in[j];
  }
}
#endif

// out[size] を 0, 1, 2, 3, ... で埋める
__global__ void kernel_iota(int* out, unsigned int size, int bias =0) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = i + bias;
  }
} 

int main() {
  using namespace std;

  const int n = 2000;
  using record = pair<int,pair<float,float>>; // key and 2 values
  vector<record> aos(n);

  mt19937 gen;
  uniform_real_distribution<float> dist;
  // 元データをつくって
  for ( int i = 0; i < n; ++i ) { aos[i] = make_pair(i, make_pair(dist(gen),dist(gen))); }
  cout << "original:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << aos[i].first << ':' 
         << aos[i].second.first << ',' << aos[i].second.second << ' ';
  } 
  cout << "...\n";

  // かき混ぜる
  shuffle(aos.begin(), aos.end(), gen);
  cout << "shuffled:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << aos[i].first << ':' 
         << aos[i].second.first << ',' << aos[i].second.second << ' ';
  } 
  cout << "...\n";

  vector<int>   h_key(n);
  vector<float> h_value1(n);
  vector<float> h_value2(n);

  // AOS から SOA にコピー
  for ( int i = 0; i < n; ++i ) {
    h_key[i]    = aos[i].first;
    h_value1[i] = aos[i].second.first;
    h_value2[i] = aos[i].second.second;
  }

  int*   d_key;
  float* d_value1;
  float* d_value2;

  float* d_value_in;
  int*   d_map;

  cudaMalloc(&d_key,      n*sizeof(int));
  cudaMalloc(&d_value1,   n*sizeof(float));
  cudaMalloc(&d_value2,   n*sizeof(float));

  cudaMalloc(&d_value_in, n*sizeof(float));
  cudaMalloc(&d_map,      n*sizeof(int));

  // keyをdeviceにコピーし、mapを初期化
  cudaMemcpy(d_key, h_key.data(), n*sizeof(int), cudaMemcpyHostToDevice);
  kernel_iota<<<(n+255)/256,256>>>(d_map, n); // 0, 1, 2, ... n-1 をセット

  // key/mapをソートする
  thrust::device_ptr<int>   p_key(d_key);
  thrust::device_ptr<int>   p_map(d_map);

  thrust::sort_by_key(p_key, p_key+n, p_map);

  // mapに基づいてvalue1を入れ替え
  cudaMemcpy(d_value_in, h_value1.data(), n*sizeof(float), cudaMemcpyHostToDevice);
#ifdef GATHER_BY_THRUST
  thrust::gather(p_map, p_map+n, 
                 thrust::device_ptr<float>(d_value_in), 
                 thrust::device_ptr<float>(d_value1));
#else
  kernel_gather<<<(n + 255)/256, 256>>>(d_map, d_value_in, d_value1, n);
#endif

  // mapに基づいてvalue2を入れ替え
  cudaMemcpy(d_value_in, h_value2.data(), n*sizeof(float), cudaMemcpyHostToDevice);
#ifdef GATHER_BY_THRUST
  thrust::gather(p_map, p_map+n, 
                 thrust::device_ptr<float>(d_value_in), 
                 thrust::device_ptr<float>(d_value2));
#else
  kernel_gather<<<(n + 255)/256, 256>>>(d_map, d_value_in, d_value2, n);
#endif

  // 結果をhostに書き戻して
  cudaMemcpy(h_key.data(),    d_key,    n*sizeof(int),   cudaMemcpyDeviceToHost);
  cudaMemcpy(h_value1.data(), d_value1, n*sizeof(float), cudaMemcpyDeviceToHost);
  cudaMemcpy(h_value2.data(), d_value2, n*sizeof(float), cudaMemcpyDeviceToHost);

  // 確認
  cout << "sorted:\n";
  for ( int i = 0; i < 10; ++i) { 
    cout << h_key[i] << ':' 
         << h_value1[i] << ',' << h_value2[i] << ' ';
  } 
  cout << "...\n";

  // あとしまつ
  cudaFree(d_key);
  cudaFree(d_value1);
  cudaFree(d_value2);
  cudaFree(d_value_in);
  cudaFree(d_map);
  cudaDeviceReset();
}



2016-06-20 07:40 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/117/

Thrust: key-value pair をソートする
2016年06月14日(火)
イマドキのプログラミング言語ではデータのまとまりを構造体で表現します:
struct record {
  int   key;
  float value;
};
record data[N]; // AOS: Array of structure
かたやGPUは飛び飛びのデータをアクセスするのが得意じゃありません(パフォーマンスが落ちちゃいます)から、データごとに配列を用意します:
  // SOA: Structure of Array
  int   key[N];
  float value[N];
このとき、ふたつの配列 key[]とvalue[]は連動させにゃならんです。
「key[]をソートする」際にはkey[]の各要素の大小に応じて要素の入れ替えを行いますが、key[i]とkey[j]を入れ替えたなら、併せてvalue[i]とvalue[j]も入れ替えないと辻褄が合わなくなりますわな。

Thrustの中におあつらえ向きのソート関数:sort_by_keyを見つけました。key[]のソートと連動してvalue[]も一緒に入れ替えてくれます。

#include <cuda_runtime.h>
#include <thrust/device_ptr.h>
#include <thrust/sort.h>

#include <iostream>
#include <vector>    // vector
#include <random>    // random
#include <utility>   // pair
#include <algorithm> // transform, shuffle

int main() {
  using namespace std;

  const int n = 2000;
  vector<pair<int,float>> aos(n);

  mt19937 gen;
  uniform_real_distribution<float> dist;
  // 元データをつくって
  for ( int i = 0; i < n; ++i ) { aos[i] = make_pair(i, dist(gen)); }
  cout << "original:\n";
  for ( int i = 0; i < 10; ++i) { cout << aos[i].first << ':' << aos[i].second << ' '; } 
  cout << "...\n";

  // かき混ぜる
  shuffle(aos.begin(), aos.end(), gen);
  cout << "shuffled:\n";
  for ( int i = 0; i < 10; ++i) { cout << aos[i].first << ':' << aos[i].second << ' '; } 
  cout << "...\n";

  vector<int>   h_key(n);
  vector<float> h_value(n);

  // AOS から SOA にコピー
  transform(aos.begin(), aos.end(), h_key.begin(),   
            [](pair<int,float> p) { return p.first;});
  transform(aos.begin(), aos.end(), h_value.begin(), 
            [](pair<int,float> p) { return p.second;});

  // key/value を deviceにコピーし
  int* d_key;
  float* d_value;
  cudaMalloc(&d_value, n*sizeof(float));
  cudaMalloc(&d_key, n*sizeof(int));

  cudaMemcpy(d_key,   h_key.data(),   n*sizeof(int),   cudaMemcpyHostToDevice);
  cudaMemcpy(d_value, h_value.data(), n*sizeof(float), cudaMemcpyHostToDevice);

  // key/valueをソートする
  thrust::device_ptr<int>   p_key(d_key);
  thrust::device_ptr<float> p_value(d_value);

  thrust::sort_by_key(p_key, p_key+n, p_value);

  // 結果をhostに書き戻して
  cudaMemcpy(h_key.data(),   d_key,   n*sizeof(int),   cudaMemcpyDeviceToHost);
  cudaMemcpy(h_value.data(), d_value, n*sizeof(float), cudaMemcpyDeviceToHost);

  // 確認
  cout << "sorted:\n";
  for ( int i = 0; i < 10; ++i) { cout << h_key[i] << ':' << h_value[i] << ' '; } 
  cout << "...\n";

  // あとしまつ
  cudaFree(d_key);
  cudaFree(d_value);
  cudaDeviceReset();
}



2016-06-14 09:40 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/116/

nvGRAPH: 最短路長(サンプルre-write)
2016年06月06日(月)
CUDA 8で新ライブラリ: nvGRAPHが追加されました。

グラフ理論に関するライブラリで、グラフすなわち
「二つの点(vertex)を結ぶ辺(edge)」の集合(topology)と各辺の重み(weights)
を食わせておくと、とある始点から各点に至る最短路長/最長路長なんかを求めてくれたりします。

で、上記topologyは疎行列で表現され、ライブラリが提供する各関数に応じてCSR/CSC-formatで用意せんなりません。

CUDA Toolkit 8RC をインストールすると nvGRAPHのサンプルコードもオマケについてくるんですけど、このサンプルがみょーにコ難しく書かれてたので、マニュアルと格闘しつつC++でre-writeしました。
/* nvGRAPH_SSSP(改)
   SSSP(Single Source Shortest Path) 
     指定した始点から各点に至る最短路長を求める
 */
#define _SCL_SECURE_NO_WARNINGS
#include <nvgraph.h>
#include <iostream>
#include <cstdlib>

void check_status(nvgraphStatus_t status) {
  if ((int)status != NVGRAPH_STATUS_SUCCESS ) {
    std::cerr << "ERROR : " << nvgraphStatusGetString(status) << std::endl;
    exit(0);
  }
}

struct edge {
  int   source;      // 始点
  int   destination; // 終点
  float weight;      // 重み(始点と終点間の長さ)
  edge(int s, int d, float w) : source(s), destination(d), weight(w) {}
  edge() {}
};

void sssp(nvgraphCSCTopology32I_t topology, const float* weights, int source_vertex, float* path_distances) {
    // nvGRAPHのハンドルを生成
    nvgraphHandle_t     handle;
    check_status(nvgraphCreate(&handle));

    // handle上にグラフを生成
    nvgraphGraphDescr_t graph;
    check_status(nvgraphCreateGraphDescr (handle, &graph));

    // - グラフにトポロジをセットし
    check_status(nvgraphSetGraphStructure(handle, graph, topology, NVGRAPH_CSC_32));

    // vertex領域をfloatでひとつ確保(ここを結果の格納領域に用いる)
    cudaDataType_t dimT = CUDA_R_32F;
    check_status(nvgraphAllocateVertexData(handle, graph, 1, &dimT));
    // edge領域をfloatでひとつ確保(ここを重みの格納領域に用いる)
    cudaDataType_t edge_dimT = CUDA_R_32F;
    check_status(nvgraphAllocateEdgeData(handle, graph, 1, &dimT));

    // - edge領域#0 に各edgeの重みをセット
    check_status(nvgraphSetEdgeData(handle, graph, (void*)weights, 0, NVGRAPH_CSC_32));

    // edge領域#0に設定された重みを使って頂点source_vertexから各点までの最短路長をvertex領域#0に求める
    check_status(nvgraphSssp(handle, graph, 0,  &source_vertex, 0));
    // vertex領域#0をpath_distancesにコピー
    check_status(nvgraphGetVertexData(handle, graph, path_distances, 0, NVGRAPH_CSC_32));

    // あとしまつ
    check_status(nvgraphDestroyGraphDescr (handle, graph));
    check_status(nvgraphDestroy (handle));
}

#include <fstream>
#include <algorithm>
#include <iterator>
#include <string>

int main() {
    using namespace std;

    const int nvertices = 6; // 頂点(vertex)の総数
    const int nedges = 10;   // 辺(edge)の総数

   // 各edgeの (始点, 終点, 重み) の列挙
    edge edges[nedges] = {
      { 0, 1, 0.500000f }, // 点0 から 点1 まで 長さ0.5 以下同文
      { 0, 2, 0.500000f },
      { 2, 0, 0.333333f },
      { 2, 1, 0.333333f },
      { 2, 4, 0.333333f },
      { 3, 4, 0.500000f },
      { 3, 5, 0.500000f },
      { 4, 3, 0.500000f },
      { 4, 5, 0.500000f },
      { 5, 3, 1.000000f }, 
    };

    if ( true ) { // オマケ: dot-formatで出力(for graphviz)
      ofstream fs("sssp.dot");
      fs <<   "digraph {\n" 
              "  node [shape=circle]\n";
      for ( edge e : edges ) {
        fs << "  " << e.source << " -> " << e.destination 
           << " [label=" << e.weight << "]\n";
      }
      fs <<   "}" << endl;
    }

// これより、上記ののCOO-format から CSC-format に変換する
// --- COO -> CSC ここから
    // 終点(destination)をmajor-key, 始点(source)をminor-keyとして比較するファンクタ
    auto comp = [](const edge& a, const edge& b) {
                    return std::make_pair(a.destination, a.source) < 
                           std::make_pair(b.destination, b.source);
                };
    // compを使ったソート
    sort(begin(edges), end(edges), comp);

    int   destination_offsets[nvertices+1]; // 圧縮された終点
    int   source_indices[nedges];           // 始点
    float weights[nedges];                  // 各edgeの重み

    // 終点列を'圧縮'して destination_offsetsに
    for ( int i = 0; i <= nvertices; ++i ) { 
      destination_offsets[i] = 
        (int)distance(begin(edges), 
                      lower_bound(begin(edges), end(edges), 
                                  edge(0, i, 0.0f), comp));
    }
    // 始点(source_indices)と重み(weights)にコピー
    transform(begin(edges), end(edges), source_indices, [](const edge& e) { return e.source;} );
    transform(begin(edges), end(edges), weights,        [](const edge& e) { return e.weight;} );
// --- COO -> CSC ここまで

    // CSR-formatを基にトポロジを構成する
    nvgraphCSCTopology32I_st topology;
    topology.nvertices           = nvertices;
    topology.nedges              = nedges;
    topology.destination_offsets = destination_offsets;
    topology.source_indices      = source_indices;

    // 計算結果をここに求める
    float path_distances[nvertices];

    // source_vertexを起点とし、各点に至る最短路長を求める
    for ( int source_vertex : { 0, 5 } ) { 
      cudaEvent_t start, stop;
      cudaEventCreate(&start); cudaEventCreate(&stop);
      cudaEventRecord(start); // startイベントを記録
      sssp(&topology, weights, source_vertex, path_distances);
      cudaEventRecord(stop);  // stopイベントを記録
      cudaEventSynchronize(stop);
      float elapsed;
      cudaEventElapsedTime(&elapsed, start, stop); // start-stop間の時間を elapsedに求める
      cout << elapsed << "[ms]\n";
      // 結果の出力
      cout << "shortest path distances from #" << source_vertex << endl;
      for ( int i = 0; i < nvertices; ++i ) {
        float val = path_distances[i];
        cout << i << " : " << ((val == FLT_MAX) ? "can't be reached" : to_string(val)) << endl;
      }
      cout << endl;
    }
    cudaDeviceReset();
}</コード>こんな結果が得られましたよ:

左側はCOO-formatから生成したdotをgraphvizでpngに変換したもの。

2016-06-06 12:25 | 記事へ | コメント(0) | トラックバック(0) |
| nvGRAPH |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/115/

cuSPARSE: cuBLASと掛け算競争
2016年05月30日(月)
cuSPARSEによる疎行列(sparse)の掛け算 とcuBLASでの密行列(dense)の掛け算、双方の処理時間を比べてみました。

1000x1000の成分をデタラメに埋めたfloat行列をふたつ用意し、その積を求めます。百万個ある成分はほとんどが0で、充填率1%(=1万)で非0値をセットしました。

GeForce GTX750, i7-2600K, PCI-gen2 での結果は:
sparse 2.0275[ms] dense 2.1668[ms]

まあ同程度の処理時間、充填率を小さくするとsparseはもっと速くなります。

ついでに GTX980, i7-6700, PCI-gen3 でもやってみたところ:
sparse 1.5534[ms] dense 0.5298[ms]

denseが速くなるのは予想どおり、GTX750と980の"格の違い"てやつなんでしょうが、対してsparseのスピードがイマイチです。

Visual Profiler で覗いたtimelineがコチラ:



cuSPARSEでの掛け算は2step(上図の青い帯)で行われます。1'st-stepの結果を基に必要なメモリを確保し、2'nd-stepで積を求めるんだけど、step間に大きなスキマがあります。このスキマはメモリの確保、つまりGPUによる計算そのものが速くなった分、メモリの確保に要する時間が相対的に大きくなったんですね。

試しにあらかじめ十分なメモリを用意して同じことやらせたら(当然ながら)このスキマが埋まり、ぐっと速くなったです。

2016-05-30 08:18 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / cuSPARSE |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/114/

CUDA Toolkit 8 RC
2016年05月28日(土)
でましたよ、CUDA Toolkit 8 RC
Release Candidate ですから、大きな不具合がなければこのまま正式版に昇格です。

Visual Studio 2015 に対応したのが Windows版の大きなトピック、待ってました! てカンジです。

それと、CUDAコンパイラ: nvcc が速くなりました。
早速インストールして 7.5 と 8RC とを比べてみたところ:


おぉ、3倍速くなってます♪

2016-05-28 12:11 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / nvcc |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/113/

cuSPARSE: 疎行列の積を求める
2016年05月23日(月)
前回のつづき、疎行列の積を求めます。

cuSPARSEにはふたつの疎行列の積を求める関数:cusparseScsrgemmが用意されています。こいつに前回使った疎行列 Aとその転置行列 tAを食わせましょうか。行列 Aは4行5列だからその転置行列 tAは5行4列、したがって両者の積は4行4列になりますな。

cusparseScsrgemm()は行列の積を求め、各成分をあらかじめ用意しておいた領域にセットしてくれるんだけど、ここで困りごとが。
"あらかじめ用意"って、どんだけ用意しとくんだ、と。

乗算の結果もまたcsr-formatで得られます。ってことは、乗算の結果にいくつの非0成分が現れるかわからんことには結果の格納領域を"あらかじめ用意"できません。

cuSPARSEでは行列の積をふたつのstepで求めます。
ひとつめのstepで積の非0成分の個数を求めてくれる関数:cusparseXcsrgemmNnzを使って必要な領域のサイズを求め、それに基づいて領域確保します。
しかるのちふたつめのstep:上述のcusparseScsrgemm()で乗算を行うってスンポーです。
#include <cuda_runtime.h>
#include <cusparse.h>
#include <algorithm>
#include <iostream>

using namespace std;

template<typename Function>
void enumerate_csr(const int* rowPtr, const int* colInd, int m, Function fn) {
  for ( int i = 0; i < m; ++i ) {
    for ( int j = rowPtr[i]; j != rowPtr[i+1]; ++j ) {
      fn(i, colInd[j], j);
    }
  }
}

int main() {
/*
  dense format
  A[m][k] = [ [ 1.0 4.0 0.0 0.0 0.0 ]
              [ 0.0 2.0 3.0 0.0 0.0 ]
              [ 5.0 0.0 0.0 7.0 8.0 ]
              [ 0.0 0.0 9.0 0.0 6.0 ] ]
  m = 4, k = 5

  B[k][n] = transpose(A), k = 5, n = 4 
  C[m][n] = A * B         m = 4. n = 4
 */

  const int m = 4;
  const int k = 5;
  const int n = 4;
  const int nnz = 9; // number of nonzero: 非0成分の総数

  float cooVal[nnz]    = { 1.0f, 4.0f, 2.0f, 3.0f, 5.0f, 7.0f, 8.0f, 9.0f, 6.0f };
  int   cooRowInd[nnz] = {   0,    0,    1,    1,    2,    2,    2,    3,    3  };
  int   cooColInd[nnz] = {   0,    1,    1,    2,    0,    3,    4,    2,    4  };

  int*  csrColInd = cooColInd; // coo と csr の ColInd は同じ

  cusparseHandle_t handle;
  cusparseCreate(&handle);
  cusparseSetPointerMode(handle, CUSPARSE_POINTER_MODE_HOST);

  cusparseMatDescr_t descr;
  cusparseCreateMatDescr(&descr); 
  cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
  cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);  

  const int nnzA = nnz;
  float* csrValA;
  int* cooRowIndA;
  int* csrRowPtrA;
  int* csrColIndA;

  cudaMallocManaged(&csrValA, nnzA*sizeof(float));
  cudaMallocManaged(&cooRowIndA, nnzA*sizeof(float));
  cudaMallocManaged(&csrRowPtrA, (m+1)*sizeof(int));
  cudaMallocManaged(&csrColIndA, nnzA*sizeof(int));

  copy(cooVal,    cooVal+nnzA,    csrValA);
  copy(cooRowInd, cooRowInd+nnzA, cooRowIndA);
  copy(csrColInd, csrColInd+nnzA, csrColIndA);

  // coo -> csr (RowInd -> RowPtr)
  cusparseXcoo2csr(handle, 
                   cooRowIndA, nnzA,
                   m, csrRowPtrA,
                   CUSPARSE_INDEX_BASE_ZERO);

  // B は A と同一
  const int nnzB  = nnzA;
  float* csrValB  = csrValA;
  int* csrRowPtrB = csrRowPtrA;
  int* csrColIndB = csrColIndA;

  int nnzC;
  float* csrValC;
  int* csrRowPtrC;
  int* csrColIndC;

  // この時点で大きさのわかっているのは 行数だけ
  cudaMallocManaged(&csrRowPtrC, (m+1)*sizeof(int));
  // 乗算の結果 非0成分の RowPtrとnnzを求めておく
  cusparseXcsrgemmNnz(handle, 
                      CUSPARSE_OPERATION_NON_TRANSPOSE,
                      CUSPARSE_OPERATION_TRANSPOSE,
                      m, n, k,
                      descr, nnzA, csrRowPtrA, csrColIndA,
                      descr, nnzB, csrRowPtrB, csrColIndB,
                      descr, csrRowPtrC, &nnzC);
  // 得られたnnz(非0成分の数)に応じてメモリ確保
  cudaMallocManaged(&csrValC, nnzC*sizeof(float));
  cudaMallocManaged(&csrColIndC, nnzC*sizeof(int));

  // C = A x transpose(B)
  cusparseScsrgemm(handle,
                   CUSPARSE_OPERATION_NON_TRANSPOSE, // A はそのまま
                   CUSPARSE_OPERATION_TRANSPOSE,     // B は転置
                   m, n, k,
                   descr, nnzA, csrValA, csrRowPtrA, csrColIndA,
                   descr, nnzB, csrValB, csrRowPtrB, csrColIndB,
                   descr,       csrValC, csrRowPtrC, csrColIndC);

  cudaDeviceSynchronize(); // 読みだす前に synchronize!!
  enumerate_csr(csrRowPtrC, csrColIndC, m, 
                [=](int r, int c, int v) {
                  cout << '[' << r << ',' << c << "] " 
                       << '\t' << csrValC[v] << endl;
                });


  // あとしまつ
  cusparseDestroyMatDescr(descr); 
  cusparseDestroy(handle);

  cudaFree(csrValA);
  cudaFree(cooRowIndA);
  cudaFree(csrRowPtrA);
  cudaFree(csrColIndA);

  cudaFree(csrValC);
  cudaFree(csrRowPtrC);
  cudaFree(csrColIndC);

}



転置行列との積なので対称行列(C[i,j] = C[j,i])ができてますです。

2016-05-23 07:59 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / cuSPARSE |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/112/

cuSPARSE: 疎行列の表現形式
2016年05月16日(月)
計算コアの数にモノを言わせ、大量のデータひとつひとつに対し単一の演算を一気に行うのがCUDAの持ち味ですから、行列演算はCUDAの得意技のひとつ。行列演算ライブラリ:cuBLASなんてのがToolkitに入っています。

ここで問題となるのがデバイスメモリの制限、CUDAが扱えるメモリはグラボに積まれた数GBのデバイスメモリなわけで、巨大な行列はデバイスメモリを圧迫します。たとえば一万行/一万列の行列だと成分数は一億、float行列ひとつで400MBを消費します。百万×百万なら4TB...これじゃグラボがパンクしますわな。

成分のほとんどが0の行列を疎行列(sparse matrix)と称します。疎行列なら非0成分のみを領域確保することで使用メモリを圧縮することができ、かなり大きな行列でもデバイスメモリに載せることができます。疎行列を対象とする行列演算ライブラリがcuSPARSEです。

で、疎行列のデバイスメモリ上での表現形式(format)を考えてみましょう、たとえばこんな4行5列の行列:
  A[4][5] = [ [ 1.0 4.0 0.0 0.0 0.0 ]
              [ 0.0 2.0 3.0 0.0 0.0 ]
              [ 5.0 0.0 0.0 7.0 8.0 ]
              [ 0.0 0.0 9.0 0.0 6.0 ] ]
成分は20個あるけど非0成分は9個です。
coo-format(Coodinate Format)で表現すると:
  const int nnz = 9; // number of nonzero: 非0成分の総数
  float cooVal[nnz]    = { 1.0f, 4.0f, 2.0f, 3.0f, 5.0f, 7.0f, 8.0f, 9.0f, 6.0f };
  int   cooRowInd[nnz] = {   0,    0,    1,    1,    2,    2,    2,    3,    3  };
  int   cooColInd[nnz] = {   0,    1,    1,    2,    0,    3,    4,    2,    4  };
cooRowInd[i]行 cooColInd[i]列 の成分が cooVal[i] であることを示しています。
ここでcooRowInd[],cooColInd[]は行index,列indexでソートしてあります。要するに行列の非0成分を第0行左から右, 第1行左から右...に列挙したやつね。

cuSPARSEの各種行列演算関数ではcoo-formatをさらに圧縮したcsr-format(Compressed Sparse Row Format)なんかが用いられます。

cuSPARSEのformat変換関数:cusparseXcoo2csr を使って coo-fromat から csr-format に変換してみましょう:
#include <cuda_runtime.h>
#include <cusparse.h>
#include <algorithm>
#include <iostream>

using namespace std;

int main() {
/*
  dense format
  A[4][5] = [ [ 1.0 4.0 0.0 0.0 0.0 ]
              [ 0.0 2.0 3.0 0.0 0.0 ]
              [ 5.0 0.0 0.0 7.0 8.0 ]
              [ 0.0 0.0 9.0 0.0 6.0 ] ]
  rows = 4, cols = 5 
 */

  const int rows = 4;
  const int cols = 5;
  const int nnz = 9; // number of nonzero: 非0成分の総数
  float cooVal[nnz]    = { 1.0f, 4.0f, 2.0f, 3.0f, 5.0f, 7.0f, 8.0f, 9.0f, 6.0f };
  int   cooRowInd[nnz] = {   0,    0,    1,    1,    2,    2,    2,    3,    3  };
  int   cooColInd[nnz] = {   0,    1,    1,    2,    0,    3,    4,    2,    4  };

  cusparseHandle_t handle;
  cusparseCreate(&handle);

  int* d_cooRowInd; cudaMallocManaged(&d_cooRowInd, nnz*sizeof(int));
  copy(cooRowInd, cooRowInd+nnz, d_cooRowInd);

  int* d_csrRowPtr; cudaMallocManaged(&d_csrRowPtr, (rows+1)*sizeof(int));

  cusparseXcoo2csr(handle, 
                   d_cooRowInd, // 行index 配列
                   nnz,         // 非0成分の数
                   rows,        // 行数
                   d_csrRowPtr, // output
                   CUSPARSE_INDEX_BASE_ZERO); // 0-base index
  cudaDeviceSynchronize();
  cout << "csrRowPtr[rows+1] = { ";
  for ( int i = 0; i < rows+1; ++i ) {
    cout << d_csrRowPtr[i] << ", ";
  }
  cout << "};" << endl;

  cudaFree(d_cooRowInd);
  cudaFree(d_csrRowPtr);
}



これがナニを意味するか、つーとですね。

「第i行の各成分は csrRowPtr[i] から csrRowPtr[i+1]に満たない範囲にある」

ことを示しています。

たとえば csrRowPtr[2] = 4, csrRowPtr[3] = 7 ですから、4から始まり7に満たない範囲:
csrRowPtr[rows+1] = { 0, 2, 4, 7, 9, };
cooVal[nnz]    = { 1.0f, 4.0f, 2.0f, 3.0f, 5.0f, 7.0f, 8.0f, 9.0f, 6.0f };
cooColInd[nnz] = {   0,    1,    1,    2,    0,    3,    4,    2,    4  };
が、第2行の非0成分です。csrRowPtr[]の大きさは(行数+1)なので、行数が非0成分の総数より小さければcoo-formatよりコンパクトですし、特定の行を参照するのにはcoo-formatより高速ってわけ。

2016-05-16 09:09 | 記事へ | コメント(0) | トラックバック(0) |
| cuSPARSE |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/111/

コンパイル時間の測定
2016年05月09日(月)
host/deviceコードのコンパイラ:nvcc はフツーのnativeコンパイラよりやんなきゃいけないことが多いんです。

まず食わせたコードをhost側とdevice側とに分離、それに併せて __global__ を呼び出してる func<<<...>>>(〜) からstubを生成してhostとdeviceを仲介します。そしてdeviceコードの生成、さらにhostコンパイラでhostコードをコンパイルののちhost/deviceそれぞれのobjを一本にマージして云々...

そんなこんなでnvccによるコンパイルは相応の時間がかかります。

NVIDIAの"なかのひと"から、マニュアルには明記されていないnvccの裏(?)オプションを教えてもらいました。

コンパイルオプションに -time <file> を指定すると、一連のコンパイル中にkickしたsub-commandと所要時間を CSV-format で出力してくれます。



出力された time.csv を excel で覗いてみると:



リリースが近い CUDA 8.0 はnvccによるコンパイルがかなり速くなるとのこと、hostコンパイラ(ここではcl.exe)が動いてる時間は変わらんでしょうから、それを除いた cudafe,cicc あたりの改善が期待できますね。

2016-05-09 08:37 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / nvcc |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/110/

cudaMemcpy2Dの使いみち
2016年04月25日(月)
__global__ void addKernel(int *c, const int *a, const int *b) {
  int i = threadIdx.x;
  c[i] = a[i] + b[i];
}
Windows(Visual C++)版 CUDAをお使いのミナサマには見飽きた感の漂うaddKernel、配列の要素ごとの加算を行うカーネルです。

こんな(よくある)お題があったとしましょうか:
struct score {
  int math;  // 数学の得点
  int engl;  // 英語の得点
  int total; // 合計点
};
...
    const int arraySize = 5;
    score exam[arraySize] = {
      { 1, 10}, { 2, 20}, { 3, 30}, { 4, 40}, { 5, 50}
    }; 
学生5人分の数学と英語のテストの得点がexam[arraySize]に納められてて、各人の合計点をexam[i].totalに求めよ、と。

このお題におなじみaddKernelを使うのは少々面倒です。struct内のメンバ:math,engl から int配列にコピーし、合計点を格納するint配列と共にaddKernelに与えて計算し、求められた合計点をexam[i].totalに書き戻すことになります。

    // allocate device memories
    int* d_a; cudaMalloc(&d_a, arraySize*sizeof(int));
    int* d_b; cudaMalloc(&d_b, arraySize*sizeof(int));
    int* d_c; cudaMalloc(&d_c, arraySize*sizeof(int));

    // allocate temp. host memories
    int* h_a = new int[arraySize];
    int* h_b = new int[arraySize];
    int* h_c = new int[arraySize];

    // copy exam[] to temp.
    for ( int i = 0; i< arraySize; ++i ) {
      h_a[i] = exam[i].math;
      h_b[i] = exam[i].engl;
    }

    // copy host to device
    cudaMemcpy(d_a, h_a, arraySize*sizeof(int), cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, arraySize*sizeof(int), cudaMemcpyHostToDevice);

    // invoke kernel
    addKernel<<<1,arraySize>>>(d_c, d_a, d_b);

    // copy device to host
    cudaMemcpy(h_c, d_c, arraySize*sizeof(int), cudaMemcpyDeviceToHost);

    // copy temp to exam
    // Add vectors in parallel.
    for ( int i = 0; i< arraySize; ++i ) {
      exam[i].total = h_c[i];
    }

    // deallocate temp 
    delete[] h_a;
    delete[] h_b;
    delete[] h_c;

    // print result
    for ( int i = 0; i< arraySize; ++i ) {
      std::cout << exam[i].math << " + "
                << exam[i].engl << " = "
                << exam[i].total << std::endl;
    }

    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);

難しくはないけどメンドくさいのに閉口しますね。device-memory d_a, d_b, d_c はCUDA使う限りやむを得ないとしても、host-memory h_a, h_b, h_c がうざったい。やりたいのは host-deviceコピーなのに host-host-deviceコピーせんならんやないですか。

初段のhost-hostコピーはstructに起因するスキマを詰める処理ですわね。こんな"スキマの空いたコピー"にcudaMemcpy2Dが使えるんですわ。

int main() {
    const int arraySize = 5;
    score exam[arraySize] = {
      { 1, 10}, { 2, 20}, { 3, 30}, { 4, 40}, { 5, 50}
    }; 

    // allocate device memories
    int* d_a; cudaMalloc(&d_a, arraySize*sizeof(int));
    int* d_b; cudaMalloc(&d_b, arraySize*sizeof(int));
    int* d_c; cudaMalloc(&d_c, arraySize*sizeof(int));

    // copy host to device
    // cudaMemcpy2D(コピー先, コピー先の幅, コピー元, コピー元の幅
    //              要素の幅, 要素数, 方向)
    cudaMemcpy2D(d_a, sizeof(int), &exam->math, sizeof(score),  
                 sizeof(int), arraySize, cudaMemcpyHostToDevice);
    cudaMemcpy2D(d_b, sizeof(int), &exam->engl, sizeof(score),  
                 sizeof(int), arraySize, cudaMemcpyHostToDevice);

    // invoke kernel
    addKernel<<<1,arraySize>>>(d_c, d_a, d_b);

    // copy device to host
    cudaMemcpy2D(&exam->total, sizeof(score), d_c, sizeof(int), 
                 sizeof(int), arraySize, cudaMemcpyDeviceToHost);

    // print result
    for ( int i = 0; i< arraySize; ++i ) {
      std::cout << exam[i].math << " + "
                << exam[i].engl << " = "
                << exam[i].total << std::endl;
    }

    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
}

2016-04-25 08:23 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/109/

線形補間(linear interpolation)
2016年04月18日(月)
等速直線運動している点があり、そのx座標(yでもzでもいいんだけどさ)の位置が 時刻 t=0.0 のとき v0, t=1.0 のとき v1 であるときに、
任意の時刻 t(0.0〜1.0) での位置 v は:

v = (1-t)*v0 + t * v1;

で求まります。線形補間(lerp)ってやつです。CUDAだと:
__device__ inline lerp(float v0, float v1, float t) {
  return (1.0f - t)*v0 + t*v1;
}
この計算は 加減算と積算をそれぞれ2回行いますし、演算のたびに結果の丸めが行われます。

CUDA Math API に fma(fused multiply-add:融合積和)てのがありましてね、fmaf(float x, float y, float z) は x*y+z を一度に計算し、丸めも一回です。上のlerpをfmaf使って書き直すと:
__device__ inline lerp(float v0, float v1, float t) {
  return fmaf(t, v1, fmaf(-t, v0, v0));
}
ほんの少し(数%)速く、誤差も小さくなるらしいですよ。

2016-04-18 09:22 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / tips |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/108/

device lambda のキャプチャ
2016年04月11日(月)
CUDA 7.5 以降、デバイス側でlambdaが使えるようになりました。

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

template<typename Function>
__global__ void kernel_do_n(Function fn, unsigned int n) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < n ) fn(i);
}

int main() {
  const unsigned int N = 10;
  kernel_do_n<<<1,N>>>([] __device__ (unsigned int i) { 
                         printf("%d ", i); 
                       }, 
                       N);
  cudaDeviceSynchronize();
  puts("");
}


kernel_do_n に与えた lambda式はデバイス側で実行されるわけですが、この device lambda はホスト側の値をキャプチャできるんやろか...やってみました。
int main() {
  const unsigned int N = 10;
  unsigned int bias = 5;
  kernel_do_n<<<1,N>>>([=] __device__ (unsigned int i) { 
                         printf("%d ", i+bias); // biasをキャプチャ
                       }, 
                       N);
  cudaDeviceSynchronize();
  puts("");
}


できてますねー。ならばこんなのはどぉよ...
int main() {
  const unsigned int N = 10;
  int* d_data;
  cudaMalloc(&d_data, N*sizeof(int));
  int val = 123;
  // d_data と val をキャプチャ
  auto fill = [=] __device__ (unsigned int i) { d_data[i] = val; };
  kernel_do_n<<<1,N>>>(fill, N);
  cudaDeviceSynchronize();

  int h_data[N];
  cudaMemcpy(h_data, d_data, N*sizeof(int), cudaMemcpyDeviceToHost);
  for ( int item : h_data ) {
    printf("%d ", item);
  }
  puts("");
  cudaFree(d_data);
}



問題なさそうです。ついでに参照キャプチャも:
int main() {
  const unsigned int N = 10;
  int* d_data;
  cudaMalloc(&d_data, N*sizeof(int));
  int val = 123;
  // d_data を(参照)キャプチャ
  auto fill = [&] __device__ (unsigned int i) { d_data[i] = 123; };
  kernel_do_n<<<1,N>>>(fill, N);
  cudaDeviceSynchronize();

  int h_data[N];
  cudaMemcpy(h_data, d_data, N*sizeof(int), cudaMemcpyDeviceToHost);
  for ( int item : h_data ) {
    printf("%d ", item);
  }
  puts("");
  cudaFree(d_data);
}


叱られました。
参照キャプチャてのはよーするにアドレスをキャプチャするわけで、ホスト側にあるモノのアドレスをデバイス側に引き込めるわけがない。

2016-04-11 08:39 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/107/

__shfl_xor (3) blockまるっとreduction
2016年04月04日(月)
CUDAは山ほどあるデータの処理が持ち味ですからね、32個のreductionじゃつまんない。
__shfl_xorを使った32個のreductionをふたつ重ねれば最大1024個のreductionができるやないの。

template<typename T>
__device__ __forceinline__ int warp_reduce(T val) {
  val += __shfl_xor(val, 16);
  val += __shfl_xor(val,  8);
  val += __shfl_xor(val,  4);
  val += __shfl_xor(val,  2);
  val += __shfl_xor(val,  1);
  return val;
}

template<typename T>
__global__ void kernel_reduce(const T* data, T* result, unsigned int size) {
  __shared__ T cache[32];
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int laneID = threadIdx.x % warpSize;
  unsigned int warpID = threadIdx.x / warpSize;
  T localSum = i < size ? data[i] : 0;
  localSum = warp_reduce(localSum);
  if ( laneID == 0 ) cache[warpID] = localSum;
  __syncthreads();
  if ( threadIdx.x < warpSize ) {
    localSum = threadIdx.x < blockDim.x / warpSize ? cache[laneID] : 0;
  }
  if ( warpID == 0 ) localSum = warp_reduce(localSum);
  if ( threadIdx.x == 0 ) atomicAdd(result, localSum);
}

__syncthreads() までは前回とほぼ同じ、ちがうところは複数のwarpに対応すべく shared memory上に配列:cache[]を用意し、n番warpが求めたredunctionの結果をcache[n]にセットしています。block内のthread数は最大1024なので、cache[]の大きさは32(=1024/32)あれば十分。

block内の全threadが__syncthreads()で待ち合わせてくれるので、ココを抜けた時点では全warpがwarp_reduce()を終え、cache[]は全部埋まっています。

あとはwarpごとのreductionの結果cache[]に対し再度waro_reduceすれば、block内の全thread数分のreductionが完了です。んでもってそいつをatomicAddで足し込んで完了、と。

2016-04-04 07:47 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / warp |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/106/

__shfl_xor (2) 32個のreduction
2016年03月30日(水)
__shfl_xor を使うと同一warp内の他のthreadと値を交換できます。コレの応用。

__device__ __forceinline__ int warp_reduce(int val) {
  val += __shfl_xor(val, 16);
  val += __shfl_xor(val,  8);
  val += __shfl_xor(val,  4);
  val += __shfl_xor(val,  2);
  val += __shfl_xor(val,  1);
  return val;
}

// data[0] + data[1] + ... + data[31] を求める
__global__ void kernel_reduce_32(const int* data, int* result) {
  int val = warp_reduce(data[threadIdx.x]);
  if ( threadIdx.x == 0 ) *result = val;
}

device関数 warp_reduce、コレを呼び出したthreadのlaneが0だとしましょう。

__shfl_xor(val, 16) はlane-16が与えたvalが返ってきますから、
val += __shfl_xor(val, 16) によって valには data[0] + data[16] が求まります。
次の __shfl_xor(val, 8) はlane-8が与えたval、すなわち data[8] + data[24] が返ってきます。なので
val += __shfl_xor(val, 8) によって valには data[0] + data[16] + data[8] + data[24] が求まります。
こうやって順次やってけば、最終的に data[0]〜data[31] の総和が求まるってスンポーです。

2016-03-30 09:52 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / warp |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/105/

__shfl_xor
2016年03月28日(月)
warp命令のおべんきょ中。

kernel_fn<<<G,B>>>(〜) すると、G*B個のthreadが動き出します。このときB(blockあたりのthread数)個のthreadはwarpと呼ばれるまとまりに小分けにされ、warpごとにStreaming Multiprocessorに投げ込まれてスケジュールされます。warpひとつが面倒見るthread数はKepler,Maxwellともに32、256threads/blockであれば 256/32 = 8個のwarpがスケジュールされるってことですな。

warpひとつに32本のlaneがあり、各laneにひとつずつ、threadが乗っかります。各threadが乗ってるlaneID(lane番号)は、一次元blockの場合 laneID = threadIdx.x % warpSize で得られます。

で、warp命令のひとつ:__shfl_xor(val, mask) が面白いふるまいをみせてくれまして、"コレを発行したthreadのlaneIDとmaskとをXORして得られたlaneID"を持つthreadが与えたvalが返ってきます。

たとえばmask=16とすると laneIDが 0番のthreadには16番, 1番には17番... が与えたvalが返ってきます。

/* shfl_xor.cu
 * nvcc -Xcompiler -wd4819 -arch sm_30 shfl_xor.cu 
 */

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

__global__ void kernel_shfl_xor() {
  int give = threadIdx.x;
  int take = __shfl_xor(give, 16);
  printf("%2d->%2d  ", give, take);
}

int main() {
  kernel_shfl_xor<<<1,64>>>();
  cudaDeviceSynchronize();
}


よーするに、同一warpに乗っかってるthread間でなら、work領域なしで値の交換ができるですよ。

2016-03-28 10:06 | 記事へ | コメント(0) | トラックバック(0) |
| warp / CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/104/

transform2D
2016年03月20日(日)
配列 X[N] と 関数 f(x) があって Y[i] = f(X[i]), i = 0,1...N-1 を求める関数 transform の二次元版を作ってみました。

配列 X[M][N] と 関数 f(x) を食わすと Y[j][i] = f(X[j][i]) を求めます。

#ifndef TRANSFORM2D_H_
#define TRANSFORM2D_H_

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

template<typename T, typename U, typename DeviceFunction>
__global__ void kernel_transform2D(const T* src, int srcPitch,
                                         U* dst, int dstPitch,
                                         int width, int height,
                                         DeviceFunction func) {
  int x = blockDim.x * blockIdx.x + threadIdx.x;
  int y = blockDim.y * blockIdx.y + threadIdx.y;
  if ( x < width && y < height ) {
    const T* pSrc = (const T*)((char*)src + srcPitch*y);
          U* pDst = (      U*)((char*)dst + dstPitch*y);
    pDst[x] = func(pSrc[x]); 
  }
}

template<typename T, typename U, typename DeviceFunction>
cudaError_t transform2D(const T* src, int srcPitch,
                              U* dst, int dstPitch,
                              int width, int height,
                              DeviceFunction func,
                              cudaStream_t stream=nullptr) {
  dim3 block(32,8);
  dim3 grid((width + 31)/32, (height + 7)/8);
  kernel_transform2D<<<grid,block>>>(src, srcPitch, dst, dstPitch, 
                                     width, height, func);
  return cudaGetLastError();
}

#endif
思いつきで書いてみたんだけども、これがなかなか使い勝手がよろしいようで、たとえばカラー画像のセピア変換とか、lambdaひとつ用意すればイッパツでできちゃうですよ。

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <opencv2/opencv.hpp>
#include <iostream>

#include "transform2D.h"

__device__ inline unsigned char clamp(float val) {
  if ( val <   0.0f ) val =   0.0f;
  if ( val > 255.0f ) val = 255.0f;
  return static_cast<unsigned char>(val);
}

using namespace std;

int main() {
  cv::VideoCapture capture(-1);
  if ( !capture.isOpened() ) {
    cerr << "can't create VideoCapture" << endl;
    return -1;
  }

  cv::String color_window  = "Color";
  cv::String sepia_window  = "Sepia";
  cv::namedWindow(color_window);
  cv::namedWindow(sepia_window);
  cv::Mat frame;

  uchar3* d_color = nullptr;
  int     color_step;

  while ( cv::waitKey(1) != 27 ) {
    // フレームをキャプチャ
    capture >> frame;
    // 元画像を表示
    cv::imshow(color_window, frame);

    IplImage image = frame;

    // device-memory を確保(only once)
    if ( d_color == nullptr ) {
      if ( image.nChannels != 3 || image.depth != 8 ) {
        cerr << "sorry, unsupported frame-format" << endl;
        break;
      }
      size_t pitch;
      cudaMallocPitch(&d_color, &pitch, image.width * sizeof(uchar3), image.height);
      color_step = static_cast<int>(pitch);
    }

    // 画像データをホストからデバイスにコピー
    cudaMemcpy2D(d_color, color_step, image.imageData, image.widthStep, 
                 image.width*sizeof(uchar3), image.height, cudaMemcpyHostToDevice);

    // カラー(BGR) を セピアに色調変換
    transform2D(d_color, color_step, d_color, color_step, image.width, image.height,
                [] __device__ (uchar3 val) -> uchar3 { 
                  float Y = val.x*0.114f + val.y*0.587f + val.z*0.299f;
                  return make_uchar3(clamp(Y-44.1915f), clamp(Y+18.8832f), clamp(Y+17.8831f));
                });

    // 処理された画像データをデバイスからホストにコピー
    cudaMemcpy2D(image.imageData, image.widthStep, d_color, color_step, 
                 image.width*sizeof(uchar3), image.height, cudaMemcpyDeviceToHost);

    // セピア画像を表示
    cv::imshow(sepia_window, frame);
  }
  cv::destroyAllWindows();

  cudaFree(d_color);
}

2016-03-20 13:57 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/103/

行列を転置する(そのに)
2016年03月14日(月)
前回 はcuBLASのヘルパ関数 cublassSetVectorを何度か呼び出して行列を転置してみたんだけど、device内で転置するカーネルを書いてみた。

ひとつは device内の行列をふたつ用意し、一方の転置を他方に作るもの。
もうひとつは行列の中で行/列を反転させるもの。

#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <assert.h>
#include "transpose.h"

template<typename T>
__global__ void kernel_transpose(int width, int height,
                                 const T* src, int ldSrc,
                                       T* dst, int ldDst) {
  unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
  if ( x < width && y < height ) {
    const T* pSrc = src + y*ldSrc + x;
          T* pDst = dst + x*ldDst + y;
    *pDst = *pSrc;
  }
}

template<typename T>
__host__ void transpose(int width, int height,
                        const T* src, int ldSrc,
                              T* dst, int ldDst,
                        cudaStream_t stream) {
  assert( width  <= ldSrc );
  assert( height <= ldDst );
  dim3 block(32,8);
  dim3 grid((width+31)/32, (height+7)/8);
  kernel_transpose<<<grid,block,0,stream>>>(width, height, src, ldSrc, dst, ldDst);
}

void transposeS(int width, int height, const float*           src, int ldSrc, float*           dst, int ldDst, cudaStream_t stream) 
  { transpose<float>(width,height,src,ldSrc,dst,ldDst,stream); }
void transposeD(int width, int height, const double*          src, int ldSrc, double*          dst, int ldDst, cudaStream_t stream)
  { transpose<double>(width,height,src,ldSrc,dst,ldDst,stream); }
void transposeC(int width, int height, const cuComplex*       src, int ldSrc, cuComplex*       dst, int ldDst, cudaStream_t stream)
  { transpose<cuComplex>(width,height,src,ldSrc,dst,ldDst,stream); }
void transposeZ(int width, int height, const cuDoubleComplex* src, int ldSrc, cuDoubleComplex* dst, int ldDst, cudaStream_t stream)
  { transpose<cuDoubleComplex>(width,height,src,ldSrc,dst,ldDst,stream); }

void transpose(int width, int height, const float*           src, int ldSrc, float*           dst, int ldDst, cudaStream_t stream) 
  { transpose<float>(width,height,src,ldSrc,dst,ldDst,stream); }
void transpose(int width, int height, const double*          src, int ldSrc, double*          dst, int ldDst, cudaStream_t stream)
  { transpose<double>(width,height,src,ldSrc,dst,ldDst,stream); }
void transpose(int width, int height, const cuComplex*       src, int ldSrc, cuComplex*       dst, int ldDst, cudaStream_t stream)
  { transpose<cuComplex>(width,height,src,ldSrc,dst,ldDst,stream); }
void transpose(int width, int height, const cuDoubleComplex* src, int ldSrc, cuDoubleComplex* dst, int ldDst, cudaStream_t stream)
  { transpose<cuDoubleComplex>(width,height,src,ldSrc,dst,ldDst,stream); }

template<typename T>
__global__ void kernel_transpose(int width, T* mtx, int ldMtx) {
  unsigned int x = blockDim.x * blockIdx.x + threadIdx.x;
  unsigned int y = blockDim.y * blockIdx.y + threadIdx.y;
  if ( x < width && y < width && x < y ) {
    T* src = mtx + y*ldMtx + x;
    T* dst = mtx + x*ldMtx + y;
    T tmp = *dst;
     *dst = *src;
     *src =  tmp;
  }
}

template<typename T>
__host__ void transpose(int width, T* mtx, int ldMtx, cudaStream_t stream) {
  assert( width  <= ldMtx );
  dim3 block(32,8);
  dim3 grid((width+31)/32, (width+7)/8);
  kernel_transpose<<<grid,block,0,stream>>>(width, mtx, ldMtx);
}

void transposeSinplace(int width, float*           mtx, int ldMtx, cudaStream_t stream)
  { transpose<float>(width, mtx, ldMtx,stream); }
void transposeDinplace(int width, double*          mtx, int ldMtx, cudaStream_t stream)
  { transpose<double>(width, mtx, ldMtx,stream); }
void transposeCinplace(int width, cuComplex*       mtx, int ldMtx, cudaStream_t stream)
  { transpose<cuComplex>(width, mtx, ldMtx,stream); }
void transposeZinplace(int width, cuDoubleComplex* mtx, int ldMtx, cudaStream_t stream)
  { transpose<cuDoubleComplex>(width, mtx, ldMtx,stream); }

void transpose(int width, float*           mtx, int ldMtx, cudaStream_t stream)
  { transpose<float>(width, mtx, ldMtx,stream); }
void transpose(int width, double*          mtx, int ldMtx, cudaStream_t stream)
  { transpose<double>(width, mtx, ldMtx,stream); }
void transpose(int width, cuComplex*       mtx, int ldMtx, cudaStream_t stream)
  { transpose<cuComplex>(width, mtx, ldMtx,stream); }
void transpose(int width, cuDoubleComplex* mtx, int ldMtx, cudaStream_t stream)
  { transpose<cuDoubleComplex>(width, mtx, ldMtx,stream); }




2016-03-14 08:03 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/102/

行列を転置する
2016年03月07日(月)
線型代数ライブラリ cuBLAS、1979年にFORTRANで実装されたBLAS(Basic Liner Algebra Subprograms)のCUDA-portです。

元がFORTRANであるがために、C/C++だとちょっと居心地の悪いとこがあります。

FORTRANでの行列すなわち二次元配列のメモリ・レイアウトはcolumn-majorであるのに対しC/C++ではrow-major、cuBLASに行列を与えるとき、行と列を転置させにゃならんのです。

たとえば行列:
| 1 2 3 |
| 4 5 6 |
だと、
float mtx[2][3] = {
  { 1, 2, 3 },
  { 4, 5, 6 }
};
じゃなくて、
float mtx[3][2] = {
  { 1, 4 },
  { 2, 5 },
  { 3, 6 }
};
のように転置させておくことになります。

この表記にはなかなか慣れないので、行列の転置関数を作ってみました。
cuBLASのヘルパ関数 cublasSetVector/cublasGetVector() がいいカンジに使えそう。

cublasSetVector(n, elemSizxe, x, incx, y, incy) は「elemSizeバイトの要素をn個、xから始まりincx個おきに読んで yから始まるincyおきに書き込」んでくれるので、incx/incyを適切に選べば 横一行 と 縦一列 の配置替えができます。
この操作を行数or列数繰り返して転置を行うのでパフォーマンスが気になりますが、host⇔deviceコピーはそんなに頻繁にやるこっちゃないでしょうから大目にみてあげますか。
# ホントに速くしたければdevice内で転置するカーネルを書くことになるかと。

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

cublasStatus_t
cublasSetMatrixTranspose(int rows, int cols, int elemSize, 
                         const void* A, int lda, void* B, int ldb) {
  cublasStatus_t status;
  if ( cols < rows ) {
    for ( int i = 0; i < cols; ++i ) {
      status = cublasSetVector(rows, elemSize, A, lda, B, 1);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + elemSize;
      B = (char*)B + ldb*elemSize; 
    }
  } else {
    for ( int i = 0; i < rows; ++i ) {
      status = cublasSetVector(cols, elemSize, A, 1, B, ldb);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + lda*elemSize;
      B = (char*)B + elemSize; 
    }
  }
  return status;
}

cublasStatus_t
cublasSetMatrixTransposeAsync(int rows, int cols, int elemSize, 
                              const void* A, int lda, void* B, int ldb, 
                              cudaStream_t stream) {
  cublasStatus_t status;
  if ( cols < rows ) {
    for ( int i = 0; i < cols; ++i ) {
      status = cublasSetVectorAsync(rows, elemSize, A, lda, B, 1, stream);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + elemSize;
      B = (char*)B + ldb*elemSize; 
    }
  } else {
    for ( int i = 0; i < rows; ++i ) {
      status = cublasSetVectorAsync(cols, elemSize, A, 1, B, ldb, stream);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + lda*elemSize;
      B = (char*)B + elemSize; 
    }
  }
  return status;
}

cublasStatus_t
cublasGetMatrixTranspose(int rows, int cols, int elemSize, 
                         const void* A, int lda, void* B, int ldb) {
  cublasStatus_t status;
  if ( cols < rows ) {
    for ( int i = 0; i < cols; ++i ) {
      status = cublasGetVector(rows, elemSize, A, lda, B, 1);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + elemSize;
      B = (char*)B + ldb*elemSize; 
    }
  } else {
    for ( int i = 0; i < rows; ++i ) {
      status = cublasGetVector(cols, elemSize, A, 1, B, ldb);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + lda*elemSize;
      B = (char*)B + elemSize; 
    }
  }
  return status;
}

cublasStatus_t
cublasGetMatrixTransposeAsync(int rows, int cols, int elemSize, 
                              const void* A, int lda, void* B, int ldb, 
                              cudaStream_t stream) {
  cublasStatus_t status;
  if ( cols < rows ) {
    for ( int i = 0; i < cols; ++i ) {
      status = cublasGetVectorAsync(rows, elemSize, A, lda, B, 1, stream);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + elemSize;
      B = (char*)B + ldb*elemSize; 
    }
  } else {
    for ( int i = 0; i < rows; ++i ) {
      status = cublasGetVectorAsync(cols, elemSize, A, 1, B, ldb, stream);
      if ( status != CUBLAS_STATUS_SUCCESS ) break;
      A = (const char*)A + lda*elemSize;
      B = (char*)B + elemSize; 
    }
  }
  return status;
}

// ----- おためし -----

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

int main() {
  cublasHandle_t handle;
  cublasCreate(&handle);

  const int ROW = 3;
  const int COL = 4;

  const int ldSrc = 128; // COL;
  const int ldDst = 64; // ROW;

  float data[ROW][ldSrc] = { { 11, 12, 13, 14 },
                             { 21, 22, 23, 24 },
                             { 31, 32, 33, 34 } };

  float* src = data[0]; // float[ROW][COL]
  float* dst;           // float[COL][ROW]
  cudaMallocManaged(&dst, COL*ldDst*sizeof(float));

  cublasSetMatrixTransposeAsync(ROW, COL, sizeof(float), src, ldSrc, dst, ldDst, nullptr);
  cudaStreamSynchronize(nullptr);

  for ( int y = 0; y < COL; ++y ) {
    for ( int x = 0; x < ROW; ++x ) {
      cout << dst[ldDst*y + x] << ' ';
    }
    cout << endl;
  }

  fill_n(data[0], ROW*COL, 0.0f);
  cublasGetMatrixTransposeAsync(COL, ROW, sizeof(float), dst, ldDst, src, ldSrc, nullptr);
  cudaStreamSynchronize(nullptr);

  for ( int y = 0; y < ROW; ++y ) {
    for ( int x = 0; x < COL; ++x ) {
      cout << data[y][x] << ' ';
    }
    cout << endl;
  }

  cudaFree(dst);
  cublasDestroy(handle);
  cudaDeviceReset();
}

2016-03-07 07:42 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / cuBLAS |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/101/

device-memoryの利用率
2016年02月29日(月)
CUDA runtime API に cudaMemGetInfo があります。
  size_t free;
  size_t total;
  cudaMemGetInfo(&free, &total);
こんなコードでdevice-memoryの空き容量/総量を取得できます。
コレを使ってですね、

1. cudaDeviceReset直後の空き容量:initial_freeを取得する
2. chunkバイトのdevice-memoryをcudaMallocで確保する
3. (2)をrep回繰り返す
4. (3)終了時点での空き容量:current_freeを取得する

そーすっと deviceから確保したメモリの総量は rep*chunkバイト、
対してdevice-memory上の使用量は (initial_free - current_free)
となり、両者の比がdevice-memoryの利用率となります。

chunk = 1, 2, ... 2^n について利用率を計測しました。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <iostream>

void measure_consumption(size_t chunk, size_t limit) {
  cudaDeviceReset();
  size_t initial_free; // 最初の容量
  size_t current_free;
  size_t total;
  cudaMemGetInfo(&initial_free, &total);
  size_t rep = 0;
  void* ptr = nullptr;
  while ( rep*chunk < limit && cudaMalloc(&ptr,chunk) == cudaSuccess ) {
    ++rep;
  }
  cudaMemGetInfo(¤t_free, &total); // 確保後の容量
  size_t actual_consume = initial_free - current_free; // 使用量
  std::cout << chunk << ", " << rep << ", " 
            << rep*chunk << ", " << (float)rep*chunk/actual_consume 
            << std::endl;
  cudaDeviceReset();
}

int main() {
  size_t chunk = 1U;
  size_t limit = 1024*1024*256; // 256MB
  for ( int i = 0; i < 24; ++i ) {
    std::cout << i << ", ";
    measure_consumption(chunk, limit);
    chunk *= 2U;
  }
}

GeForce GT720 ではこんな結果になりました(GTX750でも同様)



この結果から察するに、cudaMallocでdevice-memoryを確保するとき、どんなに小さなサイズを指定してもdevice内では512バイトを消費するようです。 ってことは、(512バイトに満たない)小さな領域をちびちび確保するより512バイト(floatなら128個分)まとめて確保し、アプリケーション側で切り分けて使った方がdevice上の無駄がないってことですな。

どうやらこの512バイトがdeviceには「キリのいい」単位のようです。
cudaMallocPitchで二次元領域を確保すると、pitcvhサイズが512の倍数になるんですよね。

2016-02-29 07:32 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/100/

shared memory のシンボルが漏れる...
2016年02月22日(月)
block内の全スレッドからアクセスできるshared memoryは、固定長/可変長の2通りの確保ができます:
  __shared__ float fcache[10]; // 固定長
  extern __shared__ double dcache[]; // 可変長

extern修飾し配列要素数を指定しなければ可変長ですな。このときshared memoryのサイズはkernel着火パラメータ<<<>>>の3つめにbyte単位で指定します。
kernel_func<<<10,10,256>>>(〜); // 256byteのshared-memory
でね、floatを引数にとるkernelを書きました:
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>

__global__ void kernel_fn(float arg) {
  extern __shared__ float cache[];
  cache[0] = arg;
  printf("kernel_fn(%f) %d:%d\n", 
         cache[0], blockIdx.x, threadIdx.x);
  /// do something
}

int main() {
  kernel_fn<<<2,3,sizeof(float) >>>(123.0f);
}
どってことなく動きます。が、これにもうひとつdouble版を追加すると:
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>

__global__ void kernel_fn(float arg) {
  extern __shared__ float cache[];
  cache[0] = arg;
  printf("kernel_fn(%f) %d:%d\n", 
         cache[0], blockIdx.x, threadIdx.x);
  /// do something
}

__global__ void kernel_fn(double arg) {
  extern __shared__ double cache[];
  cache[0] = arg;
  printf("kernel_fn(%f) %d:%d\n", 
         cache[0], blockIdx.x, threadIdx.x);
  /// do something
}

int main() {
  kernel_fn<<<2,3,sizeof(float) >>>(123.0f);
  kernel_fn<<<2,3,sizeof(double)>>>(456.00);
}
コンパイルエラーとなっちまいます。



ふたつのkernel_fnでexternしているcache[]の型が異なるぞ! だそうです。

このエラーはふたつのkernel_fnをそれぞれ別の 〜.cu に書いて コンパイル・ユニット内でユニークにしてやれば回避できます。

では二つのkernel_fnをtemplateでひとつにすると?
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>

template<typename T>
__global__ void kernel_fn(T arg) {
  extern __shared__ char cache[];
  cache[0] = arg;
  printf("kernel_fn(%f) %d:%d\n", cache[0], blockIdx.x, threadIdx.x);
  /// do something
}

int main() {
  kernel_fn<<<2,3,sizeof(float) >>>(123.0f);
  kernel_fn<<<2,3,sizeof(double)>>>(456.00);
}

するってーと先ほどのコンパイルエラーが再燃しますわね。
kernel_fn<float> と kernel_fn<double> とが異なる型の cache[] をexternしますから。

この場合コンパイル・ユニットを分けることができません、さてどうしましょ...

正解はコチラ:
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>

template<typename T>
__global__ void kernel_fn(T arg) {
  extern __shared__ char cache_[];
  T* cache = reinterpret_cast<T*>(cache_);
  cache[0] = arg;
  printf("kernel_fn(%f) %d:%d\n", 
         cache[0], blockIdx.x, threadIdx.x);
  /// do something
}

int main() {
  kernel_fn<<<2,3,sizeof(float) >>>(123.0f);
  kernel_fn<<<2,3,sizeof(double)>>>(456.00);
}

ひとまずchar配列ってことにしてコンパイラをなだめ、
しかるのち T* にキャストするってゆー、
なんとも姑息で安直な回避策ですwww

2016-02-22 08:15 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/99/

Unified Memory の舞台裏
2016年02月15日(月)
ここしばらく、Unified Memoryの挙動を観察してます。

Unified Memory(managed memmory) を確保すればその領域はhost/deviceで共有され、host-device間のコピーが不要になります。

...つっても実際にはhost-device間コピーが裏でこっそり行われてまして、プログラマが明示的にコピーせんでええ、ってことなんですが。

Visual Studio版CUDA で用意されている project:「配列の要素ごとの加算」をUnified Memoryでre-writeしました。
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <initializer_list>

cudaError_t addWithCuda(int *c, const int *a, const int *b, unsigned int size);

__global__ void addKernel(int *c, const int *a, const int *b)
{
    int i = threadIdx.x;
    c[i] = a[i] + b[i];
}

int main()
{
    const int arraySize = 5;
    int* a; cudaMallocManaged(&a, arraySize*sizeof(int));
    int* b; cudaMallocManaged(&b, arraySize*sizeof(int));
    int* c; cudaMallocManaged(&c, arraySize*sizeof(int));
    int* p;
    p = a; for ( int n : {  1,  2,  3,  4,  5 }) { *p++ = n; }
    p = b; for ( int n : { 10, 20, 30, 40, 50 }) { *p++ = n; }
    p = c; for ( int n : {  0,  0,  0,  0,  0 }) { *p++ = n; }

    // Add vectors in parallel.
    cudaError_t cudaStatus = addWithCuda(c, a, b, arraySize);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "addWithCuda failed!");
        return 1;
    }

    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 must be called before exiting in order for profiling and
    // tracing tools such as Nsight and Visual Profiler to show complete traces.
    cudaStatus = cudaDeviceReset();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceReset failed!");
        return 1;
    }

    return 0;
}

// Helper function for using CUDA to add vectors in parallel.
cudaError_t addWithCuda(int *c, const int *a, const int *b, unsigned int size)
{
    cudaError_t cudaStatus;

    // Choose which GPU to run on, change this on a multi-GPU system.
    cudaStatus = cudaSetDevice(0);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaSetDevice failed!  Do you have a CUDA-capable GPU installed?");
        goto Error;
    }

    // Launch a kernel on the GPU with one thread for each element.
    addKernel<<<1, size>>>(c, a, b);

    // Check for any errors launching the kernel
    cudaStatus = cudaGetLastError();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "addKernel launch failed: %s\n", cudaGetErrorString(cudaStatus));
        goto Error;
    }
    
    // cudaDeviceSynchronize waits for the kernel to finish, and returns
    // any errors encountered during the launch.
    cudaStatus = cudaDeviceSynchronize();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceSynchronize returned error code %d after launching addKernel!\n", cudaStatus);
        goto Error;
    }

Error:
    return cudaStatus;
}

こいつをNsight上で実行しても(裏でコッソリやってる)host-device間のコピーはtimeline上には現れません。

けどもCUDA Toolkitに同梱されているnvvp(NVIDIA Visual Profiler)なら、Unified Memoryの挙動を見せてくれます。



ふむ、kernel着火直後、配列 a,b,c がhostからdeviceへコピーされています。cは加算の結果を納める領域なのでこのコピーは不要なのですが、とにかくhostにあるものは一旦全部deviceに流し込むんですかね。

んでもって結果の読み出し、配列cを参照した時点でcをdeviceからコピーしている様子。おそらくcを参照した時点でsegment違反を検出し、それをトリガにdeviceからコピーしているんでしょう。

2016-02-15 10:11 | 記事へ | コメント(0) | トラックバック(0) |
| nvvp |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/98/

デバイスコード内の assert
2016年02月07日(日)
一連の処理が正しく進行していることの確認のため、
要所々々に assert 埋め込んだりしますわな。

CUDA ホスト側だと:
cudaError_t statusl;
...
status = cudaMemcpy(〜);
assert(status == cudaSuccess);
kernel_fun<<<grid,block>>>(〜);
assert(cudaGetLastError() == cudaSuccess);
みたいに。

このassert、デバイス側にも埋め込むことができます:
#undef NDEBUG
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <assert.h>

__global__ void kernel_test_assert(dim3 N, dim3 assertpoint) {
  dim3 global(blockDim.x*blockIdx.x+threadIdx.x, blockDim.y*blockIdx.y+threadIdx.y);
  if ( global.x < N.x && global.y < N.y ) {
    bool reached = global.x == assertpoint.x && global.y == assertpoint.y;
    printf(reached ? "*%d,%d*\t" : "(%d,%d)\t", global.x, global.y);
    assert( !reached );
  }
}

int main() {
  dim3 N(10,10);
  dim3 block(4,4);
  dim3 grid((N.x+block.x-1)/block.x, (N.y+block.y-1)/block.y);
  dim3 assertpoint(3,3); // globalIDがコレなら停止!
  cudaError_t status;

  char* h_dummy = new char[10];
  char* d_dummy; cudaMalloc(&d_dummy, 10);

  kernel_test_assert<<<grid,block>>>(N,assertpoint);
  status = cudaGetLastError();
  printf("--- launch-kernel [%s]\n", cudaGetErrorString(status));
  status = cudaMemcpy(h_dummy, d_dummy, 10, cudaMemcpyDeviceToHost);
  printf("--- cudaMemcpy    [%s]\n", cudaGetErrorString(status));
  cudaDeviceReset();
  printf("cudaDeviceReset() here\n");
  cudaMalloc(&d_dummy, 10); // リセットしちゃったのでもっぺん cudaMalloc
  status = cudaMemcpy(h_dummy, d_dummy, 10, cudaMemcpyDeviceToHost);
  printf("--- cudaMemcpy    [%s]\n", cudaGetErrorString(status));
  cudaDeviceReset();
}




ひとしきりいぢくってわかったこと:

- ひとつのスレッドでassertしても他のスレッドが道連れに停止することはないみたい。

- kernel 起動時には停止しない(非同期実行だからアタリマエ)

- assertレポートが(stderrに)出力されるのはassert発生後の同期ポイント、cudaDeviceSynchronize/cudaStreamSynchronize/cudaMemcpyなどデバイスの完了を待ち合わせるcuda-APIが呼び出されたところ。

- 一旦assertでひっかかるとその後の大抵のcuda-APIは "assertしちゃったから"って理由でハネられる。
cudaDeviceReset() でassert状態をリセットできる(devico-memoryなんかもチャラになるけど)

...まぁ、どれも納得のふるまいではありますが。


2016-02-07 09:03 | 記事へ | コメント(0) | トラックバック(0) |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/97/

managedCuda : runtime compile
2016年01月31日(日)
managedCuda、なかなかよくできてます。

前回の addKernel はあらかじめ nvcc でコンパイルして PTX を用意しておかんならんかったですが、runtime-compile使うと、kernel-codeを文字列で埋め込んでおいて実行時にコンパイルしPTXを生成できます。
/*
 * managedCuda : http://kunzmi.github.io/managedCuda/
 */
using System;
using ManagedCuda;
using ManagedCuda.NVRTC;

namespace managedCuda {
  class Program {
    static void Main(string[] args) {

      string src = 
  @"extern ""C"" __global__ 
    void addKernel(float *c, const float *a, const float *b, int size) {
      int i = blockDim.x * blockIdx.x + threadIdx.x;
      if ( i < size ) {
        c[i] = a[i] + b[i];
      }
  }";

      // 上記kernelをコンパイルしPTXを取得
      CudaRuntimeCompiler compiler = new CudaRuntimeCompiler(src, "addKernel.cu");
      compiler.Compile(null);
      byte[] ptx = compiler.GetPTX();

      // kernel(PTX中のaddKernel)をロード
      int deviceID = 0;
      CudaContext ctx = new CudaContext(deviceID);
      CudaKernel kernel = ctx.LoadKernelPTX(ptx, "addKernel");

      int N = 5;
      // grid/block 設定
      kernel.GridDimensions = (N + 255) / 256;
      kernel.BlockDimensions = 256;

      // h_A[], h_B[] をホストに確保
      float[] h_A = new float[] {  1.0f,  2.0f,  3.0f,  4.0f,  5.0f };
      float[] h_B = new float[] { 10.0f, 20.0f, 30.0f, 40.0f, 50.0f };

      // d_A[], d_B[], d_C[] をデバイスに確保
      CudaDeviceVariable<float> d_A = h_A;
      CudaDeviceVariable<float> d_B = h_B;
      CudaDeviceVariable<float> d_C = new CudaDeviceVariable<float>(N);

      // kernel着火!
      kernel.Run(d_C.DevicePointer, d_A.DevicePointer, d_B.DevicePointer,  N);

      // h_C[] を確保し 結果(d_C)をコピー
      float[] h_C = d_C;
      // コタエ合ってるかな?
      Array.ForEach(h_C, item => Console.Write("{0} ", item));
      Console.WriteLine();
    }
  }
}

2016-01-31 14:55 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / NVRTC / .NET |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/96/

managedCuda : ちょっとだけ遊んでみた
2016年01月30日(土)
CUDA-seminarでちょっとだけ紹介した namagedCuda の感触を試すべく、
Hello World 的お試しコードを書いてみました。

おなじみ配列の足し算: addKernel.cu
extern "C" __global__ 
void addKernel(float *c, const float *a, const float *b, int size) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    c[i] = a[i] + b[i];
  }
}
こいつをコンパイルし、addKernel.ptxをこしらえておきます。
> nvcc -Xcompiler -wd4819 -ptx addKernel.cu
んでもってC#コンソールアプリケーションプロジェクトを起こしてmanagedCuda.dllを参照設定に追加し、
/*
 * managedCuda : http://kunzmi.github.io/managedCuda/
 */
using System;
using ManagedCuda;

namespace managedCuda {
  class Program {
    static void Main(string[] args) {
      int N = 5;

      // kernel(addKernel.ptxの中のaddKernel)をロード
      int deviceID = 0;
      CudaContext ctx = new CudaContext(deviceID);
      CudaKernel kernel = ctx.LoadKernel("addKernel.ptx", "addKernel");

      // grid/block 設定
      kernel.GridDimensions = (N + 255) / 256;
      kernel.BlockDimensions = 256;

      // h_A[], h_B[] をホストに確保
      float[] h_A = new float[] {  1.0f,  2.0f,  3.0f,  4.0f,  5.0f };
      float[] h_B = new float[] { 10.0f, 20.0f, 30.0f, 40.0f, 50.0f };

      // d_A[], d_B[], d_C[] をデバイスに確保
      CudaDeviceVariable<float> d_A = h_A;
      CudaDeviceVariable<float> d_B = h_B;
      CudaDeviceVariable<float> d_C = new CudaDeviceVariable<float>(N);

      // kernel着火!
      kernel.Run(d_C.DevicePointer, d_A.DevicePointer, d_B.DevicePointer,  N);

      // h_C[] を確保し 結果(d_C)をコピー
      float[] h_C = d_C;
      // コタエ合ってるかな?
      Array.ForEach(h_C, item => Console.Write("{0} ", item));
      Console.WriteLine();
    }
  }
}
buildしたらば addKernel.ptx を exe-dir. にコピーして実行。


ドヤァ

2016-01-30 13:54 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/95/

nvccでのコンパイルをちょっと(かなり?)速くする
2016年01月29日(金)
CUDAセミナでは
「nvccでのコンパイルが遅いのなんとかならんか」
って質問をいただきまして。

..なんですよねー、遅いんです正直。

こんなお試しコード: hello.cu を書きました:
#include <cuda_runtime.h>
#include <stdio.h>

__global__ void kernel() {
  printf("Hello from %d\n", threadIdx.x);
}

int main() {
  kernel<<<1,5>>>();
  cudaDeviceSynchronize();
}
たったこんだけ、こいつをコンパイルします。
> mkdir tmp
> nvcc -Xcompiler -wd4819 --keep --keep-dir tmp hello.cu
コンパイルオプション:
--keep : コンパイル中に生成される中間ファイルを消さずにとっとけ
--keep-dir <dir> : 中間ファイルは<dir>に置け
を意味します。

コンパイルののちtmpのナカミを覗いてみると:



総数32本 12MB超の中間ファイルを作ってます。そらま遅いわな。

なので、この中間ファイルの置き場にSSDやramdiskとかの速いストレージをあてがうとコンパイル時間がかなり短縮されます。

上記のように --keep-dir で指定してもいいし、いちいち書くのが面倒ならば環境変数 TEMP(Windows), TMPDIR(Linux) を書き換えとくのもアレかと。

あとは〜.cuからホストコードを極力排除し __device__, __global__ と <<<うにゃうにゃ>>> なkernel着火部分だけにすることでnvccじゃないとコンパイルできない部分にダイエットするが吉でしょかね。

2016-01-29 12:30 | 記事へ | コメント(0) | トラックバック(0) |
| nvcc |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/94/

2016-01-27 CUDA seminar : aftercare
2016年01月28日(木)
先日のCUDAセミナ、たくさんのご参加ありがとうございました。

80分の長尺セッションやったんでたっぷり仕込んでおいたのですが、いざ喋ってみると"あっちゅーま"で用意してたDEMOがかなり駆け足になっちまいまして。

お見せできなかったコードがたくさんあることだし、after-careを兼ねて Visual Studio 2013 ソリューションまるごと公開します(CUDASeminarJan2016.zip)。



- CUDA.ImageConverter
  CUDAを使ったカラー画像(24bit/pix)の変換wrapper(C++/CLI)
- OpenCV.CameraCapture
  OpenCVでWeb-cameraからのフレーム・キャプチャwrapper(C++/CLI)

- OpenCV_NPP_Grayscale
  Web-cameraからのキャプチャ画像に変換施して表示(C++)
- SimpleOpenCV
  OpenCV習作:キャプチャ画像を表示するだけ(C/C++)

- OpenCV_Non_CPPCLI
  Windows Forms版 キャプチャ+変換アプリ(C++/CLI)
- WpfCameraCapture
  これがGOAL、WPF版 キャプチャ+変換アプリ(C#)

※ OpenCV 3.1 を同梱しておきました。
実行の際は OpenCV/build/x64/vc12/bin/opencv_world310.dll を exe-dir.にコピーするなりなんなり、path-reachableにしといてください。

2016-01-28 20:43 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/93/

CUDA本サンプル for Visual Studio 2013
2016年01月23日(土)

これ、販社のWeb-pageからサンプル・コードが入手できます。

各chapterごとに Linux用Makefileが入ってんですけど、Visual Studio の solution/project が含まれていません。

週末に Visual Studio 2013 solution をこしらえました。
OpenACC/MPIを使ったコードもあり、すべてのビルド/実行を確認できてはいないんですけど、WindowsでCUDAなミナサマ、お使いくださいませ。

2016-01-23 23:22 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/92/

速いGPUはメモリも速い
2016年01月18日(月)
GeForce GTX980 を手に入れました(によによ

速いっすねー、3D-mark FireStrike スコア一万を超えましたよ。
GTX980 が積んでるCUDA-coreは2048個、いつも使ってるGTX750の4倍、GT720の十倍超です。

こうなるとPCI-busを介したhost-device間のコピーにかかる時間が相対的に大きな割合を占めることになり、データ転送も速くないとせっかくの性能が薄まっちゃいます。

streamを3本立て、それぞれで host→device, kernel-call, device→host の一連の処理を非同期にやらせてみました。

まずこちらが GTX750、緑:host→device, 赤:kernel, 紫:device→host です。
非同期なのでデータ転送とkernel-callがオーバーラップしてます。



んでもってコチラが GTX980



わかります? GTX750に比べてスキマが小さい。
データ転送の上りと下りもオーバーラップしてます。

このクラスになるとコピー・エンジンをふたつ持ってて、上りと下りを同時に実行できるです。
# 聞くところによると「GTX 960以上」らしい...

2016-01-18 08:00 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/91/

atomicInc / atomicDec
2016年01月12日(火)
CUDAが提供する Lock-free アルゴリズムの中に atomicInc/atomicDec があります。いずれも特定のポインタが指すunsiged int値を +1/-1 します。

unsigned int atomicInc(unsigned int* address, unsigned int val);
unsigned int atomicDec(unsigned int* address, unsigned int val);

*addressを+1/-1してくれるワケですが、第二引数valは何なんだ、と。

第二引数valには+1/-1時の上限を指定します。
- atomicInc(address,val) の結果(*address)はvalを超えません、超えたら0に戻します。
- atomicDec(address,val) の結果(*address)は0を下回りません、底をついたらvalに戻します。

atomicInc/atomicDec の結果は 0以上・val以下 が保証される、と。

で、ですね。atomicInc(address,val)で+1 するのと atomicAdd(address,1) で+1 するの、どっちが速いか比べてみました。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <iostream>

__device__ unsigned int count;

__global__ void kernel_increment_by_atomicAdd(unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    atomicAdd(&count, 1);
  }
}

__global__ void kernel_increment_by_atomicInc(unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    atomicInc(&count, size);
  }
}

int main() {
  const unsigned int N = 10000000U;
  unsigned int tmp;

  cudaStreamQuery(nullptr);

  tmp = 0;
  cudaMemcpyToSymbol(count, &tmp, sizeof(unsigned int));
  kernel_increment_by_atomicAdd<<<(N+255)/256,256>>>(N);
  cudaMemcpyFromSymbol(&tmp, count, sizeof(unsigned int));
  std::cout << "by atomicAdd " << tmp << std::endl;

  tmp = 0;
  cudaMemcpyToSymbol(count, &tmp, sizeof(unsigned int));
  kernel_increment_by_atomicInc<<<(N+255)/256,256>>>(N);
  cudaMemcpyFromSymbol(&tmp, count, sizeof(unsigned int));
  std::cout << "by atomicInc " << tmp << std::endl;
}

Nsightでの結果がコレ:


ほんのちょびっと、atomicIncの方が速いって結果が出てますが、
1千万回やってのことですから誤差のうちじゃないかと。
同じと言っちゃっていいでしょね。



2016-01-12 11:07 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/90/

セピア変換
2016年01月07日(木)
OpenCVの準備もできたことだし、CUDA使って少しばかり気の利いたことやらせましょうかね。OpenCVのカメラ・キャプチャで得られた画像をセピア調に変換します。

cv::VideoCaptureのインスタンス:captureを作り、cv::Mat frame を用意して capture >> frame するとキャプチャされた1フレームが frame に得られます。さらに IplImage image = frame することで、image からはキャプチャされたフレームの 幅/高さ、生データの先頭アドレス等が手に入ります。これらに基づいたフレーム・データをCUDA側すなわちdevice-memoryに流し込み、CUDAで処理しようって魂胆です。

カラー画像を一旦Grayscale(モノクロ)に変換し、続いて色調変換を施してセピア調に仕立てます。

Web-cameraからキャプチャしたカラー画像は1pixelあたりBlue/Green/Red各8bit(unsigned char)の24bit(3byte)、この3つを1:6:3の比で混ぜるとイイカンジのGrayscaleになります。この処理はCUDAライブラリのひとつ、画像処理が得意なNPPにやらせます。

得られたGrayscale画像はpixelあたり8bitの輝度情報です。YUV→BGR変換において、Y:輝度とし U,Vを定数とすれば色調を変えることができます。U=0.05, V=-0.1 とするといわゆるセピア調になるらしい。この変換は小さなkernel書いてCUDAにやらせます。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <nppi.h>

__device__ inline unsigned char clamp(float val) {
  if ( val <   0.0f ) val =   0.0f;
  if ( val > 255.0f ) val = 255.0f;
  return static_cast<unsigned char>(val);
}

// 輝度Y に U, V で補正することで Grayscaleに色調変換を施す
__global__ void kernel_gray2sepia(const unsigned char* src, int src_step, 
                                        unsigned char* dst, int dst_step, 
                                  NppiSize roi, float U, float V) {
  int ix = blockDim.x * blockIdx.x + threadIdx.x;
  int iy = blockDim.y * blockIdx.y + threadIdx.y;
  if ( ix < roi.width && iy < roi.height ) {
    const float Y  = (src + src_step*iy)[ix];
    ((uchar3*)(dst + dst_step*iy))[ix] = make_uchar3(
      clamp(Y                 + 1.7330f*V*255),
      clamp(Y - 0.3444f*U*255 - 0.7114f*V*255 ),
      clamp(Y + 1.4026f*U*255));
  }
}

__host__ void launch_gray2sepia(const unsigned char* src, int src_step, 
                                   unsigned char* dst, int dst_step, 
                                   NppiSize roi, float U, float V) {
  dim3 block(32,8);
  dim3 grid((roi.width + block.x -1)/block.x, (roi.height + block.y -1)/block.y);
  kernel_gray2sepia<<<grid,block>>>(src, src_step, dst, dst_step, roi, U, V);
}

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace std;

int main() {
  cv::VideoCapture capture(-1);
  if ( !capture.isOpened() ) {
    cerr << "can't create VideoCapture" << endl;
    return -1;
  }

  cv::String color_window  = "Color";
  cv::String sepia_window  = "Sepia";
  cv::namedWindow(color_window);
  cv::namedWindow(sepia_window);
  cv::Mat frame;

  unsigned char* d_color = nullptr;
  unsigned char* d_gray  = nullptr;
  int      color_step;
  int      gray_step;
  NppiSize roi;  // Region Of Interest: 処理対象となる領域

  while ( cv::waitKey(1) != 27 ) {
    // フレームをキャプチャ
    capture >> frame;
    // 元画像を表示
    cv::imshow(color_window, frame);

    IplImage image = frame;

    // device-memory を確保
    if ( d_color == nullptr ) {
      if ( image.nChannels != 3 || image.depth != 8 ) {
        cerr << "sorry, unsupported frame-format" << endl;
        break;
      }
      roi.width = image.width;
      roi.height = image.height;
      size_t pitch;
      cudaMallocPitch(&d_color, &pitch, roi.width * sizeof(uchar3),        roi.height);
      color_step = static_cast<int>(pitch);
      cudaMallocPitch(&d_gray,  &pitch, roi.width * sizeof(unsigned char), roi.height);
      gray_step = static_cast<int>(pitch);

      cout << "width = "      << image.width     << "\theight = "   << image.height
           << "\tstride = "   << image.widthStep << "\nchannels = " << image.nChannels
           << "\tdepth = "    << image.depth     << endl;
      cerr << "*** Press ESC to exit. ***" << endl;
    }

    // 画像データをホストからデバイスにコピー
    cudaMemcpy2D(d_color, color_step, image.imageData, image.widthStep, 
                 roi.width*sizeof(uchar3), roi.height, cudaMemcpyHostToDevice);

    // カラー(BGR) を Grayscaleに変換
    Npp32f coeffs[3] = { 0.114f, 0.587f, 0.299f };
    nppiColorToGray_8u_C3C1R(d_color, color_step, d_gray, gray_step, roi, coeffs);
    // Grayscale を セピアに色調変換
    launch_gray2sepia(d_gray, gray_step, d_color, color_step, roi, 0.05f, -0.1f);

    // 処理された画像データをデバイスからホストにコピー
    cudaMemcpy2D(image.imageData, image.widthStep, d_color, color_step, 
                 roi.width*sizeof(uchar3), roi.height, cudaMemcpyDeviceToHost);

    // セピア画像を表示
    cv::imshow(sepia_window, frame);
  }
  cv::destroyAllWindows();

  cudaFree(d_color);
  cudaFree(d_gray);
}


できあがったのがコチラ:



変換にかかる時間は GT720 でもミリ秒を切る程度なんでカクカクしたりはしませんです。

2016-01-07 09:31 | 記事へ | コメント(0) | トラックバック(0) |
| OpenCV / NPP / CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/89/

Dynamic Parallelism プロパティ設定
2016年01月05日(火)
Dynamic parallelism に必要なオプションおぼえがき。

1. nvcc : Relocatable Device Code を生成させること。



2. nvcc : ComputeCapability 3.5 以上にすること。



3. link : cudadevrt.lib をリンクすること。


2016-01-05 09:10 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/88/

Dynamic Parallelism
2016年01月04日(月)
Compute Capability 3.5 以降、Dynamic Parallelismとかいう機能がCUDAに追加されています。

平たく言えば kernelを入れ子にできる(kernelからkernelが呼べる)ってことみたいで、くだらんサンプル書いてお正月を過ごしてました。

float配列 data[size] の総和を求めるkernelを用意しました:
__global__ void kernel_accumulate(float* result, const float* data, unsigned int size) {
  __shared__ float sum;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( threadIdx.x == 0 ) sum = 0.0f;
  syncthreads();
  if ( i < size ) {
    atomicAdd(&sum, data[i]);
  }
  syncthreads();
  if ( threadIdx.x == 0 ) atomicAdd(result, sum);
}
これを使って2次元配列:data[row][col]の総和を求めよう、と。
kernel_accumulateが(1次元)配列の総和を求めてくれるんだから、こいつを data[0], data[1], ..., data[row-1] に対して呼べばいいんぢゃね?
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

#include <iostream>
#include <vector>
#include <algorithm>
#include <random>
#include <numeric>
#include <cassert>

using namespace std;

__global__ void kernel_accumulate(float* result, const float* data, unsigned int size) {
  __shared__ float sum;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( threadIdx.x == 0 ) sum = 0.0f;
  syncthreads();
  if ( i < size ) {
    atomicAdd(&sum, data[i]);
  }
  syncthreads();
  if ( threadIdx.x == 0 ) atomicAdd(result, sum);
}

__global__ void kernel_accumulate2D(float* result, const float* data, unsigned int row, unsigned int col, size_t pitch) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < row ) {
    /* Dynamic Palallelism : kernel が kernel を呼ぶ! */
    kernel_accumulate<<<(col+255)/256, 256>>>(result, (const float*)((const char*)data + i*pitch), col);
  }
}

__host__ float accumulate2D(const float* data, unsigned int row, unsigned int col) {
  float* d_result;
  cudaMalloc(&d_result, sizeof(float));
  cudaMemset(d_result, 0, sizeof(float));

  float* d_data;
  size_t pitch;
  cudaMallocPitch(&d_data, &pitch, col*sizeof(float), row);
  cudaMemcpy2D(d_data, pitch, data, col*sizeof(float), col*sizeof(float), row, cudaMemcpyHostToDevice);

  kernel_accumulate2D<<<(row+255)/256, 256>>>(d_result, d_data, row, col, pitch);
  cudaDeviceSynchronize();

  float result;
  cudaMemcpy(&result, d_result, sizeof(float), cudaMemcpyDeviceToHost);

  cudaFree(d_result);
  cudaFree(d_data);
  return result;
}

int main() {
  const unsigned int R = 100;
  const unsigned int C = 100;
  vector<float> data(R*C);
  iota(begin(data), end(data), 1.0f); // 1, 2, ... N
  shuffle(begin(data), end(data), mt19937());
  
  float result;

  result = accumulate2D(data.data(), R, C);
  cout << result << endl;
}

2016-01-04 15:21 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/87/

新春にセミナーやりますよ
2015年12月29日(火)
「C#, CUDA, Industry色濃いめ」ってお題でセッション・スピーカを依頼されました。

CUDA アンバサダ の看板揚げてるからには引き受けますともさ。
NVIDIA主催、無償のCUDAセミナです、
お誘いあわせのうえご参加くださいませ。


MaxwellとJava、C#のためのCUDA
2016-01-27 13:00〜 ベルサール神田

2015-12-29 18:06 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / C++/CLI / C# |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/86/

OpenCV : 画像処理の前準備
2015年12月28日(月)
えー、今回CUDAの出番はありません ゴメンネー

「画像処理はCUDAの華」ですからね、動画の処理をやらせてみよかとWeb-cameraを買ってきました。1280x720pixのお安いやつ。

こいつを使って「フレームをキャプチャし、CUDAで処理して表示」を繰り返せばパラパラ漫画方式の動画処理ができるやろ、と。

で、まずはカメラ画像1フレームをキャプチャせにゃならんので、使い勝手の良さそうなライブラリを探し回り、画像処理ライブラリのド定番: OpenCV に決定。先日リリースされたばっかりの version 3.1 Windows版には vc12,vc14(Visual Studio 2013,2015)のbuild済ライブラリが入ってて、インストールするだけで使えます。

さっそくお味見、Webカメラからキャプチャしたフレームを表示するだけのアプリを書きました。

#if 1 /* C++ */

#include <opencv2/opencv.hpp>

int main() {
  cv::VideoCapture capture(0);
  if ( !capture.isOpened() ) {
    return -1;
  }

  cv::namedWindow("Camera Capture");
  while ( cv::waitKey(1) != 27 ) {
    cv::Mat frame;
    capture >> frame;
    cv::imshow("Camera Capture", frame);
  }
}

#else /* C */

#include <opencv/highgui.h>

int main() {
  CvCapture* capture = cvCreateCameraCapture(0);
  if ( capture == nullptr ) {
    return -1;
  }

  cvNamedWindow("Camera Capture",1);
  while ( cvWaitKey(1) != 27 ) {
    IplImage* frame = cvQueryFrame(capture);
    cvShowImage("Camera Capture", frame);
  }
  cvReleaseCapture(&capture);
}
#endif
コンパイル/実行すると:



ほんの十数行書くだけであっさりできちゃいました。

2015-12-28 21:19 | 記事へ | コメント(0) | トラックバック(0) |
| OpenCV |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/85/

nvprof の(おちゃらかな)使い方
2015年12月21日(月)
CUDA ToolkitをインストールするとVisual profiler: NsightがVisual Studio にビルトインされますけど、もうひとつ、コマンドラインから使うprofiler: nvprof もお手軽でオススメ。

コード中で #include <cuda_profiler_api.h> し、profileを停止する箇所に cudaProfilerStop(); を書き加えておきます。
んでもって nvprof --print-gpu-summary なんちゃら.exe すると host-device間のcudaMemcpy および kernel-call の回数と所要時間とを出力してくれます。

"Thrust の stream 実行" で書いたおためしコードを nvprofに食わせてみました。



kernel-callと同程度の時間を host-device 間のcopyに消費しているのが見えますね。

テキストの羅列なのでとっつきは良くないのですが、時間のかかってる順に出力されるのでボトルネックをざっくり探すのに重宝します。

# 以降に続く呪文みたいな文字列が、呼び出してるkernel関数です。
# thrust使うとこんな関数呼び出すのねー...


2015-12-21 08:44 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/84/

thrust の stream 実行
2015年12月14日(月)
Thrustが提供する関数がGPU上で実行されるとき、
デフォルトでは第0 stream(default-stream) が実行streamとなります。

が、多くのthrust関数は第1引数に「どこ(どのstream)で実行するか」を明示的に指定できます。

以前実装したサンプル:三角関数表 をthrustにportし、timelineを眺めてみます。

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

#include <thrust/device_ptr.h>
#include <thrust/transform.h>
#include <thrust/system/cuda/execution_policy.h>

#include <iostream>
using namespace std;

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

__device__ __forceinline__ float deg2rad(float angle) { return angle * 3.1416f / 180.0f; }

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

  // streams
  cudaStream_t ssin, scos, stan;
  for ( cudaStream_t* pstream : { &ssin, &scos, &stan })
    cudaStreamCreate(pstream);

  // host memories
  float* angle = nullptr; // in
  float* sin   = nullptr; // out
  float* cos   = nullptr; // out
  float* tan   = nullptr; // out
  for ( float** phost : { &angle, &sin, &cos, &tan })
    cudaMallocHost(phost, BN);

  // device memories
  float* dangle;
  float* dsin;
  float* dcos;
  float* dtan;
  for ( float** pdevice : { &dangle, &dsin, &dcos, &dtan })
    cudaMalloc(pdevice, BN);

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

  thrust::device_ptr<float> pangle(dangle);
  thrust::device_ptr<float> psin(dsin);
  thrust::device_ptr<float> pcos(dcos);
  thrust::device_ptr<float> ptan(dtan);

  cudaEvent_t trigger;
  cudaEventCreate(&trigger);

  // copy angle[] from host to device
  cudaMemcpyAsync(dangle, angle, BN, cudaMemcpyHostToDevice);
  cudaEventRecord(trigger);

  cudaStreamWaitEvent(ssin, trigger, 0);
  thrust::transform(thrust::cuda::par.on(ssin), pangle, pangle+N, psin, 
                    [] __device__ (float angle) { return sinf(deg2rad(angle)); });
  cudaStreamWaitEvent(scos, trigger, 0);
  thrust::transform(thrust::cuda::par.on(scos), pangle, pangle+N, pcos, 
                    [] __device__ (float angle) { return cosf(deg2rad(angle)); });
  cudaStreamWaitEvent(stan, trigger, 0);
  thrust::transform(thrust::cuda::par.on(stan), pangle, pangle+N, ptan,
                    [] __device__ (float angle) { return tanf(deg2rad(angle)); });

  cudaMemcpyAsync(sin, dsin, BN, cudaMemcpyDeviceToHost, ssin);
  cudaMemcpyAsync(cos, dcos, BN, cudaMemcpyDeviceToHost, scos);
  cudaMemcpyAsync(tan, dtan, BN, cudaMemcpyDeviceToHost, stan);

  cudaDeviceSynchronize();

  for ( unsigned int i = 0; i < 5; ++i ) {
    cout << angle[i] << " :"
         << " sin=" << sin[i] 
         << " cos=" << cos[i] 
         << " tan=" << tan[i] 
         << endl; 
  }
  cout << "..." << endl;

  // あとしまつ
  cudaEventDestroy(trigger);
  for ( cudaStream_t stream : { ssin, scos, stan })
    cudaStreamDestroy(stream);
  for ( float* phost : { angle, sin, cos, tan })
    cudaFreeHost(phost);
  for ( float* pdevice : { dangle, dsin, dcos, dtan })
    cudaFree(pdevice);
}

default-streamで実行したのがコレ:


対して3つのstreamで実行したときのtimelineは:


キチキチに詰まって速くなってます。

2015-12-14 08:20 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/83/

thrust::find_if
2015年12月10日(木)
STL-likeにCUDAが使えるライブラリ:Thrustにもfind_ifが用意されています。自作find_if とスピード競争してみました。

#include <thrust/device_ptr.h> // thrust::device_ptr
#include <thrust/find.h>       // thrust::find_if

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

template<typename T,
         typename Predicate>
__global__ void kernel_find_if(unsigned int* result, 
                        const T* data, 
                        const Predicate& predicate, 
                        unsigned int size) {
  __shared__ unsigned int cache[1024];
  unsigned int thread_id = threadIdx.x;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  // data[i]が条件を満たせば i, さもなくば size を cache にセットする
  cache[thread_id] = (i < size && predicate(data[i])) ? i : size;
  syncthreads();
  // cahce[] 内で最も小さな値を cache[0] に追い込む
  for ( unsigned int stride = blockDim.x/2; stride > 0; stride >>= 1) {
    if ( thread_id < stride && cache[thread_id] > cache[thread_id + stride] ) {
      cache[thread_id] = cache[thread_id + stride];
    }
    syncthreads();
  }
  // cache[0] と *result の小さい方を *result へ
  if ( thread_id == 0 ) atomicMin(result, cache[0]);
}

/* こっからおためし */

#include <iostream>
#include <vector>
#include <numeric>
#include <random>
#include <algorithm>
using namespace std;

template<typename T> inline void set(T* item, T val) 
  { cudaMemcpy(item, &val, sizeof(T), cudaMemcpyHostToDevice); }

template<typename T> inline void get(const T* item, T& val) 
  { cudaMemcpy(&val, item, sizeof(T), cudaMemcpyDeviceToHost); }

__device__ float d_val; // 関数スコープ内には置けない。残念。

int main() {

  const unsigned int N = 2000000;
  vector<float> data(N);

  cudaEvent_t start, stop;
  cudaEventCreate(&start);
  cudaEventCreate(&stop);

  float* d_data;
  cudaMalloc(&d_data, N*sizeof(float));

  unsigned int* d_result;
  cudaMalloc(&d_result, sizeof(unsigned int));

  mt19937 rnd;
  cout << "size, kernel_find_if, thrust::find_if" << endl;  
  for ( unsigned int size = 10; size < N; size *= 1.2 ) {
    cout << size << ',';
    iota(begin(data), begin(data)+size, -1.0f);
    shuffle(begin(data), begin(data)+size, rnd);
    cudaMemcpy(d_data, data.data(), size*sizeof(float), cudaMemcpyHostToDevice);

    unsigned int index;
    float elapsed;

    dim3 block = 512U; // 2のベキであること! 絶対!
    dim3 grid  = (size+block.x -1U)/block.x;

    set(d_result, size);
    { float t = 0.0f; cudaMemcpyToSymbol(d_val,&t,sizeof(float)); }
    cudaEventRecord(start);
    kernel_find_if<<<grid,block>>>(d_result, d_data,  
                                  [] __device__ (float x) { return x < d_val;},
                                  size);
    cudaEventRecord(stop); cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsed, start, stop);
    get(d_result, index);
    cerr << "kernel_find_if\t" << index << ':' << data[index] << '\t'
         << elapsed << "[ms]\n";
    cout << elapsed << ',';

    thrust::device_ptr<float> p_data(d_data);
    thrust::device_ptr<float> p_result;

    cudaEventRecord(start);
    p_result = thrust::find_if(p_data, p_data+size, 
                               [] __device__ (float x) { return x < d_val;});
    cudaEventRecord(stop); cudaEventSynchronize(stop);
    cudaEventElapsedTime(&elapsed, start, stop);
    cerr << "thrust::find_if\t" << thrust::distance(p_data, p_result) 
         << ':' << *p_result << '\t'
         << elapsed << "[ms]\n";
    cout << elapsed << endl;
  }

  cudaEventDestroy(start);
  cudaEventDestroy(stop);

  cudaFree(d_data);
  cudaFree(d_result);
}

実行結果をCSVにre-directし、Excelでグラフ描いたのがコチラ。



自作find_ifの処理時間は要素数に比例し、一直線に並んでいます。
一方 thrst::find_if は要素数の小さいときに暴れてますがこちらもおおむねリニアに増加、要素数が10万を超えたあたりで追い越します。

ちくしょー、負けた。

2015-12-10 08:10 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/82/

find_if を作ってみた
2015年12月07日(月)
reduce関数の応用、data[size] の中から特定の条件を満たす最初の要素の添え字、つまり目的の要素は何番目にあるかを返す関数に仕立ててみました。

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

template<typename T,
         typename Predicate>
__global__ void kernel_find_if(unsigned int* result, 
                        const T* data, 
                        const Predicate& predicate, 
                        unsigned int size) {
  __shared__ unsigned int cache[1024];
  unsigned int thread_id = threadIdx.x;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  // data[i]が条件を満たせば i, さもなくば size を cache にセットする
  cache[thread_id] = (i < size && predicate(data[i])) ? i : size;
  syncthreads();
  // cahce[] 内で最も小さな値を cache[0] に追い込む
  for ( unsigned int stride = blockDim.x/2; stride > 0; stride >>= 1) {
    if ( thread_id < stride && cache[thread_id] > cache[thread_id + stride] ) {
      cache[thread_id] = cache[thread_id + stride];
    }
    syncthreads();
  }
  // cache[0] と *result の小さい方を *result へ
  if ( thread_id == 0 ) atomicMin(result, cache[0]);
}

/* こっからおためし */

#include <iostream>
#include <vector>
#include <numeric>
#include <random>
#include <algorithm>
using namespace std;

template<typename T> inline void set(T* item, T val) 
  { cudaMemcpy(item, &val, sizeof(T), cudaMemcpyHostToDevice); }

template<typename T> inline void get(const T* item, T& val) 
  { cudaMemcpy(&val, item, sizeof(T), cudaMemcpyDeviceToHost); }

__device__ float d_val; // 関数スコープ内には置けない。残念。

int main() {

  const unsigned int N = 100000;
  vector<float> data(N);
  iota(begin(data), end(data), 0.0f);
  data[0] = -1.0f;
  data[1] = -1.0f;
  data[2] = (float)(N*2);
  data[3] = (float)(N*2);

  float* d_data;
  cudaMalloc(&d_data, N*sizeof(float));

  unsigned int* d_result;
  cudaMalloc(&d_result, sizeof(unsigned int));
  
  for ( int i = 0; i < 3; ++i ) {
    shuffle(begin(data), end(data), mt19937());
    cudaMemcpy(d_data, data.data(), N*sizeof(float), cudaMemcpyHostToDevice);

    float* result;

    dim3 block = 512U; // 2のベキであること! 絶対!
    dim3 grid  = (N+block.x -1U)/block.x;

    result = find_if(data.data(), data.data()+N, [](float n) { return n < 0.0f;});
    cout << "std::find_if\t" << distance(data.data(), result) << ':' << *result << endl;

    unsigned int index;

    set(d_result, N);
    { float t = 0.0f; cudaMemcpyToSymbol(d_val,&t,sizeof(float)); }
    kernel_find_if<<<grid,block>>>(d_result, d_data,  
                                  [] __device__ (float x) { return x < d_val;},
                                  N);
    get(d_result, index);
    cout << "kernel_find_if\t" << index << ':' << data[index] << endl;

    result = find_if(data.data(), data.data()+N, [=](float n) { return n > (float)N;});
    cout << "std::find_if\t\t" << distance(data.data(), result) << ':' << *result << endl;

    set(d_result, N);
    { float t = (float)(N+1); cudaMemcpyToSymbol(d_val,&t,sizeof(float)); }
    kernel_find_if<<<grid,block>>>(d_result, d_data,  
                                  [] __device__ (float x) { return x > d_val;},
                                  N);
    get(d_result, index);
    cout << "kernel_find_if\t\t" << index << ':' << data[index] << endl;

    cout << endl;

  }

  cudaFree(d_data);
  cudaFree(d_result);
}


ただねー...これがもんのすごく速いかって言われると正直自信ありません。data[0]から順に検索すれば見つかったところで処理終了なんですが、CUDAでは頭から順に処理されるって保証がないので全部舐めるっきゃないんですよ。

2015-12-07 08:28 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/81/

イロイロと速く求める
2015年11月30日(月)
前回、shared memoryを使って総和を求めました。やってたのは:

data[0] + data[1] + data[2] + ... + data[size-1]

をいくつかの部分に切り分け、各部分の和を足しこむ操作です。
これで総和が求まるのは、加算 a + b に以下のルールが成り立つからです。

■結合則 : ( a + b ) + c = a + ( b + c )
コレが成り立つおかげで、テキトーな部分に切り分けて求めることができます。

■交換則 : a + b = b + a
コレが成り立つおかげで、どんな順序で加算しても同じ結果が得られます。

■単位元の存在 : 0 + a = a
加算には単位元0が存在します。sizeがスレッド数の倍数でないとき、この単位元0を埋めることで実装を楽にしています。

結合則,交換則,単位元の存在,この3つが成り立つ演算なら、加算じゃなくても同じアルゴリズムで求めることができます。たとえば全データの積:

data[0] * data[1] * data[2] * ... * data[size-1]

あるいは全データの最大値を求める、これはすっごく小さな値を単位元にすりゃえぇです。

んなわけで総和を求めた


こいつの黄色いトコを差し替えられるようtemplateで書き換えました。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

template<typename Selection>
__device__ float atomicExch(float* address, float val, const Selection& selection) {
  int* address_as_int = (int*)address;
  int  old = *address_as_int;
  int  assumed;
  do {
    assumed = old;
    old = atomicCAS(address_as_int, assumed, 
                    __float_as_int(selection(val, __int_as_float(assumed))));
  } while ( assumed != old );
  return __int_as_float(old);
} 

template<typename Selection>
__device__ double atomicExch(double* address, double val, const Selection& selection ) {
  unsigned long long int* address_as_ull = (unsigned long long int*)address; 
  unsigned long long int old = *address_as_ull;
  unsigned long long int assumed; 
  do { 
    assumed = old; 
    old = atomicCAS(address_as_ull, assumed, 
                    __double_as_longlong(selection(val ,__longlong_as_double(assumed))));
  } while (assumed != old); 
  return __longlong_as_double(old); 
}

template<typename T,
         typename Function>
__global__ void kernel_reduce(T* result, 
                        const T* data, 
                              T  oob, 
                        const Function& function, 
                        unsigned int size) {
  __shared__ T cache[1024];
  unsigned int thread_id = threadIdx.x;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  cache[thread_id] = (i < size) ? data[i] : oob;
  syncthreads();
  for ( unsigned int stride = blockDim.x/2; stride > 0; stride >>= 1) {
    if ( thread_id < stride ) {
      cache[thread_id] = function(cache[thread_id], cache[thread_id + stride]);
    }
    syncthreads();
  }
  if ( thread_id == 0 ) atomicExch(result, cache[0], function);
}

#include <iostream>
#include <vector>
#include <numeric>
#include <random>
#include <algorithm>
#include <limits>
using namespace std;

template<typename T>
inline void clear(T* item) {
  cudaMemset(item, 0, sizeof(T)); 
}

template<typename T>
inline void set(T* item, T val) {
   cudaMemcpy(item, &val, sizeof(T), cudaMemcpyHostToDevice); 
}

template<typename T> 
inline T get(const T* item) { 
   T t; 
   cudaMemcpy(&t, item, sizeof(T), cudaMemcpyDeviceToHost); 
   return t;
}

int main() {

  const unsigned int N = 50000;
  vector<float> data(N);

  float* d_data;
  cudaMalloc(&d_data, N*sizeof(float));

  float* d_result;
  cudaMalloc(&d_result, sizeof(float));
  
  mt19937 gen;
  uniform_real_distribution<float> dist;
  generate(begin(data), end(data), [&]() { return dist(gen);});

  for ( int i = 0; i < 1; ++i ) {
    cudaMemcpy(d_data, data.data(), N*sizeof(float), cudaMemcpyHostToDevice);

    float result;

    result = accumulate(begin(data), end(data), 0.0f);
    cout << "std::accumulate\t\t\t" << result << endl;

    dim3 block = 512U; // 2のベキであること! 絶対!
    dim3 grid  = (N+block.x -1U)/block.x;

    clear(d_result);
    kernel_reduce<<<grid,block>>>(d_result, d_data, 0.0f, 
                                  [] __device__ (float x, float y) 
                                      { return x + y; },
                                  N);
    result = get(d_result);
    cout << "kernel_reduce(accumulate)\t" << result << endl;
    cout << endl;

    result = *std::min_element(begin(data), end(data));
    cout << "*std::min_element\t" << result << endl;

    float oob;
    oob = numeric_limits<float>::max();
    set(d_result, oob);
    kernel_reduce<<<grid,block>>>(d_result, d_data, oob, 
                                  [] __device__ (float x, float y) 
                                      { return x < y ? x : y; },
                                  N);
    result = get(d_result);
    cout << "kernel_reduce(min)\t" << result << endl;

    result = *std::max_element(begin(data), end(data));
    cout << "*std::max_element\t" << result << endl;

    oob = numeric_limits<float>::min();
    set(d_result, oob);
    kernel_reduce<<<grid,block>>>(d_result, d_data, oob, 
                                  [] __device__ (float x, float y) 
                                      { return x < y ? y : x; },
                                  N);
    result = get(d_result);
    cout << "kernel_reduce(max)\t" << result << endl;

    cout << endl;
    
  }

  cudaFree(d_data);
  cudaFree(d_result);
}


kernel_reduce一本で総和/最小値/最大値を求めています。
lambda式大活躍です。


2015-11-30 08:56 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/80/

総和を(速く)求める
2015年11月22日(日)
data[N] の総和: data[0] + data[1] + ... + data[N-1]
を求めるのはデータ処理のキホンのひとつ。これを実装するおハナシをatomic関数に書きました。

atomicAddを使って総和を求めるkernelを(チョー手抜きで)書くと:
__global__ void kernel_accumulate(float* result, const float *data, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    atomicAdd(result, data[i]);
  }
}
コレね、結果が正しいのはいいとして、実はそんなに速くはないんです。atomicAddは「データ競合で誤った加算をやっちまったら、もっぺんやり直す」を繰り返す戦術です。そうすることでmutexやcritical-sectionなどによる排他制御なしにデータ競合を回避します。

なので、atomicAddが各スレッドで"まばら"に呼ばれるならいいんだけど、多数のスレッドが一斉にatomicAddすると加算のやり直しが頻発して遅くなっちゃうですよ。

改良版がコチラ:
__global__ void kernel_accumulate_shared(float* result, const float *data, unsigned int size) {
  extern __shared__ float cache[];
  unsigned int thread_id = threadIdx.x;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  cache[thread_id] = (i < size) ? data[i] : 0.0f;
  syncthreads();
  for ( unsigned int stride = blockDim.x/2; stride > 0; stride >>= 1) {
    if ( thread_id < stride ) {
      cache[thread_id] += cache[thread_id + stride];
    }
    syncthreads();
  }
  if ( thread_id == 0 ) atomicAdd(result, cache[0]);
}


■ extern __shared__ cache[];
ブロック内の全スレッドが共有する cache[] を用意します。cache[]の要素数はスレッド数と同じになるようkernel呼び出し時に設定します。

■ cache[thread_id] = (i < size) ? data[i] : 0.0f;
ブロック内のスレッド数(blockDim.x)を8とすると、スレッド#0がcache[0], #1がcache[1],...#7がcache[7]をセットするので、これでcache[]は全部埋まります。これを全部加えればブロック内の総和になりますな。

■次に続くfor-loop、ちょいとばかりややこしい。
for-loopの一発目、ブロック数が8なので stride = 4 です。
次の if ( thread_id < stride ) に飛び込むのは thread_id = 0,1,2,3、それぞれのスレッドが:

スレッド#0: cache[0] = cache[0] + cache[0+4]
スレッド#1: cache[1] = cache[1] + cache[1+4]
スレッド#2: cache[2] = cache[2] + cache[2+4]
スレッド#3: cache[3] = cache[3] + cache[3+4]
ってわけで、cache[0〜3] に総和が集約されます。

loopの2発目ではさらに半分、stride = 2 ですから
スレッド#0: cache[0] = cache[0] + cache[0+2]
スレッド#1: cache[1] = cache[1] + cache[1+2]
よーするに cache[0〜1] に総和が集約されます。

これを stride > 0 の間繰り返せば最終的には cache[0]に総和が求まります。
各loop内でcache[]の添え字がカブらないのでデータ競合は起こりません。

■ if ( thread_id == 0 ) atomicAdd(result, cache[0])
最後にスレッド#0がブロックを代表してresultに足しこんで完了、と。

Nsightに処理時間を測ってもらいました。
上段が"手抜き版"、下段が"改良版"です。

ざっくり2倍速いすね♪

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

__global__ void kernel_accumulate(float* result, const float *data, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    atomicAdd(result, data[i]);
  }
}

__global__ void kernel_accumulate_shared(float* result, const float *data, unsigned int size) {
  extern __shared__ float cache[];
  unsigned int thread_id = threadIdx.x;
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  cache[thread_id] = (i < size) ? data[i] : 0.0f;
  syncthreads();
  for ( unsigned int stride = blockDim.x/2; stride > 0; stride >>= 1) {
    if ( thread_id < stride ) {
      cache[thread_id] += cache[thread_id + stride];
    }
    syncthreads();
  }
  if ( thread_id == 0 ) atomicAdd(result, cache[0]);
}

#include <iostream>
#include <vector>
#include <numeric>
using namespace std;

template<typename T>
inline void clear(T* item) {
  cudaMemset(item, 0, sizeof(T)); 
}

template<typename T> 
inline T get(const T* item) { 
   T t; 
   cudaMemcpy(&t, item, sizeof(T), cudaMemcpyDeviceToHost); 
   return t;
}

int main() {

  const unsigned int N = 1000000;
  vector<float> data(N);

  float* d_data;
  cudaMalloc(&d_data, N*sizeof(float));

  float* d_result;
  cudaMalloc(&d_result, sizeof(float));
  
  for ( int i = 0; i < 1; ++i ) {
    fill(begin(data), end(data), 1.0f/N);
    cudaMemcpy(d_data, data.data(), N*sizeof(float), cudaMemcpyHostToDevice);

    float result;

    result = accumulate(begin(data), end(data), 0.0f);
    cout << "std::accumulate\t\t\t" << result << endl;

    dim3 block = 512U; // 2のベキであること! 絶対!
    dim3 grid  = (N+block.x -1U)/block.x;
    unsigned int shared_size = block.x * sizeof(float);

    clear(d_result);
    kernel_accumulate<<<grid,block>>>(d_result, d_data, N);
    result = get(d_result);
    cout << "kernel_accumulate\t\t" << result << endl;

    clear(d_result);
    kernel_accumulate_shared<<<grid,block,shared_size>>>(d_result, d_data, N);
    result = get(d_result);
    cout << "kernel_accumulate_shared\t" << result << endl;
    cout << endl;
  }

  cudaFree(d_data);
  cudaFree(d_result);
}

2015-11-22 20:55 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/79/

Thrust: gather/scatter
2015年11月17日(火)
Thrustにちょいと気の利いた関数:gather/scatterを見つけました。

いずれも配列 in[] から out[] に値をコピーする関数なんですが、
頭から順にではなく、飛び飛びのindexでコピーしてくれます。

たとえば in[] の 3,2,5 番目をout[0..2]にコピーするなら
int map[] = { 3, 2, 5 } をgather(かき集め)に与え、

in[0..2] を out[] の 3,2,5番目にコピーするなら
上記のmapを scatter(バラ撒き)に与えます。

void gather_scatter() {
  const unsigned int N = 5;

  // 5つの複素数が、それぞれの実部と虚部に分かれて格納されてる
  float real_imag[N*2] = { 
    0.0f, 0.2f, 0.4f, 0.6f, 0.8f, // 実部
    0.1f, 0.3f, 0.5f, 0.7f, 0.9f  // 虚部
  };
  thrust::device_vector<float> d_real_imag(real_imag, real_imag+N*2);

  // かき集め/バラ撒き用index列
  int map[N*2] = { 0, 5, 1, 6, 2, 7, 3, 8, 4, 9};
  thrust::device_vector<int> d_map(map, map+N*2);

  thrust::device_vector<float> d_out(N*2);

  // gather:かき集めて実部,虚部の順に並び替える
  thrust::gather(d_map.begin(), d_map.end(), d_real_imag.begin(), d_out.begin());
  // できたかな?
  for ( int i = 0; i < N; ++i ) {
    std::cout << d_out[i*2+0] << '+' << d_out[i*2+1] << 'i' << std::endl;
  }

  // 一旦0クリア
  thrust::fill(d_real_imag.begin(), d_real_imag.end(), 0.0f);

  // scatter:バラ撒いて元に戻す
  thrust::scatter(d_out.begin(), d_out.end(), d_map.begin(), d_real_imag.begin());
  // できたかな?
  for ( float item : d_real_imag ) {
    std::cout << item << ' ';
  }
  std::cout << std::endl;
}


で、応用。R行C列の行列を転置し、C行R列に変換します。

template<typename Iter>
void print(Iter iter, int row, int column) {
  for ( int y = 0; y < row; ++y ) {
    for ( int x = 0; x < column; ++x ) {
      std::cout << *iter++ << ' ';
    }
    std::cout << std::endl;
  }
}

void transpose() {
  using namespace std;
  const int R = 3;
  const int C = 4;

  // 対応表を作る
  vector<int> tmap(R*C);
  iota(begin(tmap), end(tmap), 0);
  transform(begin(tmap),end(tmap),begin(tmap),
            [=](int i) { return i%R*C + i/R;});
  thrust::device_vector<int> d_map = tmap;

  // R行C列の行列
  vector<int> in = { 10, 11, 12, 13,
                     14, 15, 16, 17,
                     18, 19, 20, 21 };
  thrust::device_vector<int> d_in = in;
  print(begin(d_in), R, C);

  thrust::device_vector<int> d_out(C*R);

  //             このmapに基づいて          ここから     ここに 
  thrust::gather(begin(d_map), end(d_map), begin(d_in), begin(d_out));
  cout << "transposed:" << endl;
  print(begin(d_out), C, R);
}


2015-11-17 11:41 | 記事へ | コメント(0) | トラックバック(0) |
| Thrust |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/78/

CUDA on WPF(C#) つけたし
2015年11月09日(月)
WPFでマンデルブロ では、CUDAによるマンデルブロ集合の計算をひっきりなしに呼び出しては描画させてます。

このひっきりなしを実現する方法。

たとえば"描け!"ボタンを押したときに描画するのなら、"描け!"ボタンのClickイベントを捕まえて:
private void 描けボタン_Click(object sender, System.Windows.RoutedEventArgs e) {
  描く();
}


ってやるだけなんやけど、ひっきりなしに描画するには描く()を絶え間なく呼び出し続けにゃあきません。
そこで使えるのが System.Threading による非同期実行。非同期実行なので"ウラ"で勝手に仕事してくれます。
    private void おしごと() {
      描く();
      Dispatcher.BeginInvoke(new System.Action(おしごと), 
                             System.Windows.Threading.DispatcherPriority.Background,
                             null);
    }

これでおっけー、おしごと()を一回呼べば、おしごと()は自分自身を何度も(際限なく)呼び出します。

おしごと()は非同期にウラで実行されるので、スライダーぐりぐりいぢってもスピードが落ちません♪

# ...コードセット、欲しい方いらっしゃいますぅ?

2015-11-09 11:39 | 記事へ | コメント(0) | トラックバック(0) |
| .NET / C# |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/77/

CUDA on WPF(C#) おまけ
2015年11月05日(木)
で、CUDA.BipmapOpとまったく同じカラクリでマンデルブロ集合を描いてみた。



1フレームにつき26万個の複素数に対してがりがり演算するので相当にタフなお仕事なんですが、しょぼちんグラボ GT720でも133FPS 出てるくせしてCPUは余裕しゃくしゃく♪

2015-11-05 16:57 | 記事へ | コメント(0) | トラックバック(1) |
| CUDA / WPF |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/76/

CUDA on WPF(C#) そのに
2015年11月04日(水)
さてと。

CUDA.BitmapOp は.NET(=managed)から呼べにゃならんし、実装にはCUDA(=native)を利用するのでC++/CLIでCLRクラスライブラリをこしらえることになります。

CUDA.BitmapOpにはコンストラクタ/デストラクタ/ファイナライザとfillメソッドを実装します。

bitmapの大きさ(幅とか高さとか)はコンストラクタで設定するとして、fill()に引き渡すのは GBRA(4-byte)で現した塗りつぶしパターンと塗りつぶし対象のbitmap、要はbitmapをGBRAで埋め尽くせばいい。
// CUDA.BipmapOp.Assembly.h

#pragma once

namespace CUDA {

  public ref class BitmapOp {
  private:
    int width_;   // 幅(pix)
    int height_;  // 高さ(pix)
    int stride_;  // 1-line分のサイズ(byte)
    unsigned char* dev_bmp_; // device-memory

  public:
    BitmapOp(int width, int height, int stride);
   ~BitmapOp();
   !BitmapOp();

    void fill(cli::array<unsigned char>^ bgra, cli::array<unsigned char>^ bitmap);
  };

}

kernel関数ごりごり書いても大した手間ではないけども、CUDAライブラリ:NPP(NVIDIA Performance Primitives)におあつらえの関数を見つけたので。こいつ使って楽をします。

CUDA.BitmapOpの内部(private member)には device-memoryを保持し、コンストラクタで確保/デストラクタ(ファイナライザ)で解放しましょうか。

fillメソッドは至極単純:NPPの塗りつぶし関数を呼び、その結果(塗りつぶされたdevice-mmory)を引数でもらったbitmapにコピーするだけ。

// CUDA.BitmapOp.Assembly.cpp

#include "stdafx.h"

#include <cuda_runtime.h>
#include <nppi.h>

#include "CUDA.BipmapOp.Assembly.h"

namespace CUDA {

  BitmapOp::BitmapOp(int width, int height, int stride) 
    : width_(width), height_(height), stride_(stride) {
    unsigned char* ptr;
    cudaMalloc(&ptr, stride_ * height_);
    dev_bmp_ = ptr;
  }

 BitmapOp::~BitmapOp() { this->!BitmapOp(); }
 BitmapOp::!BitmapOp() { cudaFree(dev_bmp_); }

 void BitmapOp::fill(cli::array<unsigned char>^ bgra, 
                     cli::array<unsigned char>^ bitmap) {
    // ピン留めを忘れずに!
    pin_ptr<unsigned char> pin_bgra(&bgra[0]);
    pin_ptr<unsigned char> pin_bitmap(&bitmap[0]);
    NppiSize roi = { width_, height_ };
    nppiSet_8u_C4R(pin_bgra, dev_bmp_, stride_, roi);
    cudaMemcpy(pin_bitmap, dev_bmp_, stride_*height_, cudaMemcpyDeviceToHost);
 }
}

2015-11-04 13:23 | 記事へ | コメント(0) | トラックバック(1) |
| CUDA / C++/CLI / NPP |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/75/

CUDA on WPF(C#) そのいち
WPFアプリケーションでCUDAを使うココロミ。

C#・WPFアプリケーション プロジェクトのひな型を用意し、Windowの中にScrollViewer、さらにその中にImage(1024x1024pix.)。それとSliderを3つ配置します。


Imageと3つのSliderにはそれぞれ imgPanel,sldR,sldG,sldB と命名しました。
Sliderの値に応じてImageを塗りつぶすアプリに仕立てますよ。



WindowのLoadedイベント,およびSliderのValueChangedイベントにハンドラを追加し、XAMLから自動生成されたコード:MainWindow.xaml.csに手を加えてこんなの↓ができました。

using System.Windows;
using System.Windows.Media;
using System.Windows.Media.Imaging;

namespace WpfColorPanel {
  /// <summary>
  /// MainWindow.xaml の相互作用ロジック
  /// </summary>
  public partial class MainWindow : Window {

    private byte[] bitmap_; // 生のbitmap
    private int    width_;  // 幅(pix)
    private int    height_; // 高さ(pix)
    private int    stride_; // 1-line分のサイズ(byte)

    CUDA.BitmapOp bitmapop_; // bitmap操作を行うクラス

    public MainWindow() {
      InitializeComponent();
      // メンバの初期化
      width_ = (int)imgPanel.Width;
      height_ = (int)imgPanel.Height;
      stride_ = width_ * 4;
      bitmap_ = new byte[stride_*height_];
      bitmapop_ = new CUDA.BitmapOp(width_, height_, stride_);
    }

    private void fillBmp() {
      // Sliderの各値から BGRA(4-byte) を作る(Aは0固定)
      byte[] bgra = new byte[] { (byte)sldB.Value, (byte)sldG.Value, (byte)sldR.Value, 0 };
      // bitmapをbgraで塗りつぶす
      bitmapop_.fill(bgra, bitmap_);
      // bitmapからBitmapSourceを生成し、Image.Sourceにセットする
      imgPanel.Source = 
        BitmapSource.Create(width_, height_, 96, 96, 
                            PixelFormats.Bgr32, null, 
                            bitmap_, stride_);
    }

    private void Window_Loaded(object sender, RoutedEventArgs e) {
      fillBmp();
    }

    private void sld_ValueChanged(object sender, RoutedPropertyChangedEventArgs<double> e) {
      fillBmp();
    }

  }

}

ここまではどってことないですわね、あとはCUDA.BitmapOpを用意すればよし。

...長くなるので一旦ここまで。

2015-11-04 10:34 | 記事へ | コメント(0) | トラックバック(1) |
| C# / WPF |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/74/

Unified Memory 対応 vector
2015年10月26日(月)
- device-memory を確保して
- host から device に cudaMemcpy して
- kernel 呼び出して
- device から host に cudaMemcpy する

ってゆー CUDAルーチンワーク(?)が煩わしく感じることがしばしばあります。既存のコードをCUDAに置き換えてるときとか、パフォーマンスは二の次でともかくちゃんと動くか確認したいだけなんだってば! みたいなね。

Unified Memory(managed memory)を確保/解放の対象とした最小限のallocatorを用意し、Unified Memory版 vectorを作ってみたり。

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

void safe_call(cudaError_t err, const char* file, int line) {
  if ( err != cudaSuccess ) {
    std::string msg = file;
    msg += "(" + std::to_string(line) + "): ";
    msg += cudaGetErrorString(err);
    throw std::runtime_error(msg);
  }
}

#define CUDA_SAFE_CALL(x) safe_call(x, __FILE__, __LINE__)

#pragma region minimum allocator for Unified Memory
// C++11 std::allocator_traits required

template<typename T>
struct managed_allocator {
  using value_type = T;

  managed_allocator() {}
  template<typename U> managed_allocator(managed_allocator<U>&) {}

  T* allocate(std::size_t n) {
    T* ptr = nullptr;
    CUDA_SAFE_CALL(cudaMallocManaged(&ptr, n*sizeof(T)));
    return ptr;
  }

  void deallocate(T* ptr, std::size_t) {
    cudaFree(ptr);
  }

};

template<typename T, typename U>
bool operator==(const managed_allocator<T>&, const managed_allocator<U>&) { return true; }

template<typename T, typename U>
bool operator!=(const managed_allocator<T>&, const managed_allocator<U>&) { return false; }
#pragma endregion

#pragma region templatized std::vector for Unified Memory
#include <vector>
template<typename T> using managed_vector = std::vector<T,managed_allocator<T>>;
#pragma endregion


コレひとつあれば、CUDAプログラミングがぐっとお手軽になりますです。

template<typename T, typename Function>
__global__ void kernel_transform(T* out, const T* in, Function fun, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    out[i] = fun(in[i]);
  }
}

int main() {
  using namespace std;
  const unsigned int N = 1000;
  managed_vector<int> in(N);
  managed_vector<int> out(N);

  iota(begin(in), end(in), 0);

  unsigned int block = 256;
  unsigned int grid  = (N + block -1U) / block;
  kernel_transform<<<grid,block>>>(out.data(), in.data(), 
     [] __device__ (unsigned int x) -> int { return x+1; } , N);
  cudaDeviceSynchronize(); // コレ忘れずに!
//*/
  for ( int item : out ) {
    cout << item << ' ';
  }
//*/
  cout << endl;
}


2015-10-26 12:06 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / C++ |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/73/

次へ