BLOG引越のおしらせ
2016年10月09日(日)
このBlogをホストしてくれてるプロバイダさんが今年いっぱいを目途にBlogサービスを終了するってんですよ。
だもんでそろそろお引越しをはじめます。

以降のエントリは「はてな」に用意した東方算程譚に移動し、ココに残したエントリは少しずつ新居に運び込みますね。

2016-10-09 20:54 | 記事へ | コメント(0) | トラックバック(0) |
| misc. |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/133/

GTC Japan 2016 after-care
2016年10月07日(金)
GTC Japan 2016 へのご来場、ありがとうございました。

おかげさまで開発者トラックは盛況で、僕のセッションも300名は入ろうかという会場の席を埋めることができました。

ってことでセッションの最後にちょろっとお見せした
「世間の狭さをGPUが実証する」のソースコードをここに置きました。

キモとなるnvGRAPHの他、cuRAND, cuSPARSE, Thrust さらに自前のkernelと、様々なCUDAの寄せ集めで組み上げてます。

2016-10-07 22:17 | 記事へ | コメント(0) | トラックバック(0) |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/132/

GTC Japan 2016
2016年10月03日(月)
今週水曜(2016-10-05)、GTC Japan 2016 が開催されます。

僕もCUDA関連セッションをひとつ引き受けました(2016-1014)。
CUDA 8.0 で新たにお目見えした nvGRAPH をお題に25分間ヨタこきます。

いやー、セッション資料が出来上がった途端に CUDA 8.0 が正式リリースされ、RC版との相違がないか検証と修正で土日がつぶれました。

けどもおかげでドキュメントには書かれていないちょっとしたtipsを披露できそうです。

当日は会場内を回遊しているでしょうが、Hangout界隈が存在確率高メになろうかと。
見かけたら声かけてくださいねー♪

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

Streamの優先順位(そのに)
2016年09月26日(月)
んではStreamにつけられた優先順位の違いによってGPUのふるまいがどう変わるのか観察してみましょ。

二本のStreamそれぞれに、"Host→Device, kernel-call, Device→Host" の一連の命令を5組ほど投げ込みむスレッドを起こします。コード長めでゴメン。
#include <cuda_runtime.h>
#include <device_launch_parameters.h>

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

#include <iostream>
#include <thread>

struct params {
  cudaStream_t stream;
  unsigned int size;
  float* angle;
  float* sine;
  float* d_angle;
  float* d_sine;
};

void setup(params& p, unsigned int size, int priority) {
  p.size = size;
  cudaMallocHost(&p.angle, size*sizeof(float));
  cudaMallocHost(&p.sine,  size*sizeof(float));
  cudaMalloc(&p.d_angle, size*sizeof(float));
  cudaMalloc(&p.d_sine, size*sizeof(float));
  cudaStreamCreateWithPriority(&p.stream, cudaStreamDefault, priority);
}

void teardown(params& p) {
  cudaFree(p.angle);
  cudaFree(p.sine );
  cudaFree(p.d_angle);
  cudaFree(p.d_sine );
  cudaStreamDestroy(p.stream);
}

int main() {
  using namespace std;
  int low, high;
  cudaDeviceGetStreamPriorityRange(&low, &high);
//  high = low;

  const unsigned int N = 100000;
  params p0, p1;
  setup(p0, N, low);
  setup(p1, N, high);

  auto thread_fun = [](params* p) {
    for ( int i = 0; i < 5; ++i ) {
      cudaMemcpyAsync(p->d_angle, p->angle, p->size, cudaMemcpyHostToDevice, p->stream);
      kernel_sine<<<(p->size+255)/256,256,0,p->stream>>>(p->angle, p->sine, p->size);
      cudaMemcpyAsync(p->sine, p->d_sine, p->size, cudaMemcpyHostToDevice, p->stream);
    }
  };

  thread thr0(thread_fun, &p0);
  thread thr1(thread_fun, &p1);

  thr0.join(); thr1.join();

  teardown(p0);
  teardown(p1);
  cudaDeviceReset();

}
Streamの優先順位を同じにしたときのtimelineは:



まぁ、そうなるわな。

対して low と high に設定すると...



ほれ、途中でひっくり返ってますでしょ。優先順位の高いStreamから先に処理されたんですね。

その結果、後から着火したスレッドの方が先に完了するですね。

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

Streamの優先順位
2016年09月20日(火)
CUDAはキホン非同期処理。cudaMemcpyAsyncやkernel-call等はその処理が終わらぬうちにCPUに制御が戻ってきます...ってことはつまりGPUが一連の命令を溜めこむキュー(待ち行列)を持ってるってことで、それがStreamですわね。

で、Streamは(defaultの他に)複数個作ることができ、その中からお好みのStreamに命令を投げ込むことができます。

GPU(てゆーかGegaThread EngineとCopy Engine)は各Streamの先頭から命令を取り出して処理するんですが、どのStreamから命令を取り出すかはGPUの勝手です。

とはいえ、優先的に処理してくれるStreamが欲しいこともあるでしょう。

cudaStreamCreateWithPriority(cudaStream_t* pStream, unsigned int flags, int priority)
の第三引数:priority で優先順位を付けることができます。

このpriority、値が小さいほど優先順位が高く、priorityに指定できる値の範囲は
int low, high;
cudaDeviceGetPriorityRange(&low, &high);
で得られます。

lowより大きい/highより小さい値を指定するとそれぞれlow/highで頭打ち/底打ちになるみたいですよ。

※ 手持ちの GeForce GT720 では low = 0, high = -1 の二段階でした。

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

Unified Memory: nvvpで覗いてみると...
2016年09月12日(月)
先日(2016-09-10)のわんくま同盟横浜勉強会では
「CUDAめっちゃ速い...けどめっちゃ速く動かすにはそれなりに工夫しなきゃいけないのよ」
ってなセッションをやりました。

セッション終了後、"Unified Memory" のからくりに関して質問があり、よさげなデモの用意がなくてゴメンナサイしちゃったので、ここでちょろっと紹介。

float配列 angle[N], sine[N], cosine[N] を用意して、angle[i]のsin/cosを求めてsine[i]/cosine[i]に押し込むコードです:
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <math_functions.h>

__global__ void kernel_sincos(const float* angle, float* sine, float* cosine, unsigned int size) {
  int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) {
    sincosf(angle[i], sine+i, cosine+i);
  }
}

#include <iostream> // cout, endl
#include <numeric>  // iota

int main() {
  using namespace std;

  const unsigned int N = 1000000;

  float* angle;  // 角度(input)
  float* sine;   // 正弦(output)
  float* cosine; // 余弦(output)

  // 領域確保
  cudaMallocManaged(&angle,  N*sizeof(float));
  cudaMallocManaged(&sine,   N*sizeof(float));
  cudaMallocManaged(&cosine, N*sizeof(float));

  // angle[] を 0.0f, 1.0f, ... で埋める
  iota(angle, angle+N, 0.0f);

  // kernel-call
  kernel_sincos<<<(N+255)/256,256>>>(angle, sine, cosine, N);
  cudaDeviceSynchronize(); // お忘れなく!

  // こたえあわせ: sin^2 + cos^2 = 1
  float sum = 0.0f;
  for ( int i = 0; i < N; ++i ) {
    sum += sine[i]*sine[i] + cosine[i]*cosine[i];
  }
  cout << sum / N << endl;

  // あとしまつ
  cudaFree(angle);
  cudaFree(sine);
  cudaFree(cosine);
  cudaDeviceReset();
}
Unified Memory使ってますから cudaMemcpyによるhost-device間のメモリ・コピーがコード中には現れません。
そのかわり、kernel呼び出しのあと、その結果を読み出すに先立って cudaDeviceSynchronize()でdevice側の処理が完了するまで待たにゃならんです。

NVVPでこいつを動かしtimelineを覗いてみると:



kernelに着火した時点でhost→device, こたえあわせ中にdevice→host方向にメモリ・コピーが(裏でコッソリ)行われているのがバレバレです。

さらにこのtimelineをぐぐーっとZoom-inすると:



メモリ・コピーがコマ切れになってます。このコマ切れの間にはわずかながらのスキマがあるので、パフォーマンスをキッチキチに研ぎ澄ましたいのであればUnified Memoryは"向いてない"ことになります(ほんのちょびっと、ですけどね)。

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

NPP: 最小値/最大値
2016年09月05日(月)
NPPには前回のnppsSetみたいな"地味に便利"なAPIがいろいろと用意されてるみたいで、配列中から最小値と最大値を同時に求めてくれるnppsMinMaxもそのひとつ。

コレを使うには、あらかじめ作業領域をデバイス上に確保しておかにゃならんです。作業領域の大きさはnppsMinMaxGetBufferSizeで得られます。
#include <cuda_runtime.h>
#include <npps.h>

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

using namespace std;

int main() {
  const unsigned int N = 1000000;

  // 範囲 -10.0f〜10.0f の一様乱数で埋める
  vector<float> host(N);
  {
    mt19937 gen;
    uniform_real_distribution<float> dist(-10.0f, 10.0f);
    generate_n(begin(host), N, [&]() { return dist(gen); });
    auto expected = minmax_element(begin(host), end(host));
    cout << "expected: min:" << *expected.first << '\t' 
         <<           "max:" << *expected.second << endl;
  }

  // デバイス領域を確保してコピー
  float* device; 
  cudaMalloc(&device, N*sizeof(float));
  cudaMemcpy(device, host.data(), N*sizeof(float), cudaMemcpyHostToDevice);

  // 処理に必要なバッファを確保
  Npp8u* d_buffer;
  {
    int buf_size;
    nppsMinMaxGetBufferSize_32f(N, &buf_size);
    cudaMalloc(&d_buffer, buf_size);
  }

  // 結果を受け取る領域
  float* d_result; 
  cudaMalloc(&d_result, 2*sizeof(float));
  float result[2];

  // 最大値/最小値を求め
  nppsMinMax_32f(device, N, d_result+0, d_result+1, d_buffer);

  // 結果をホストにコピー
  cudaMemcpy(&result, d_result, 2*sizeof(float), cudaMemcpyDeviceToHost);

  // 結果を出力
  cout << "actual:   min:" << result[0] << '\t' 
       <<           "max:" << result[1] << endl;

  // あとしまつ
  cudaFree(device);
  cudaFree(d_buffer);
  cudaFree(d_result);
  cudaDeviceReset();
}
最小値/最大値だけでなく、それぞれの位置(index)もついでに教えてくれるnppsMinMaxIndxてーのもありますよ。
  ...
  // 処理に必要なバッファを確保
  Npp8u* d_buffer;
  {
    int buf_size;
    nppsMinMaxIndxGetBufferSize_32f(N, &buf_size);
    cudaMalloc(&d_buffer, buf_size);
  }

  // 結果を受け取る領域
  float* d_result; 
  cudaMalloc(&d_result, 2*sizeof(float));
  float result[2];
  int* d_index;
  cudaMalloc(&d_index, 2*sizeof(int));
  int index[2];

  // 最大値/最小値を求め
  nppsMinMaxIndx_32f(device, N, d_result+0, d_index+0, d_result+1, d_index+1, d_buffer);

  // 結果をホストにコピー
  cudaMemcpy(&result, d_result, 2*sizeof(float), cudaMemcpyDeviceToHost);
  cudaMemcpy(&index,  d_index,  2*sizeof(int),   cudaMemcpyDeviceToHost);

  // 結果を出力
  cout << "actual:   min:[" << index[0] << "] " << result[0] << '\t' 
       <<           "max:[" << index[1] << "] " << result[1] << endl;

  // あとしまつ
  cudaFree(device);
  cudaFree(d_buffer);
  cudaFree(d_result);
  cudaFree(d_index);
  cudaDeviceReset();
}


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

NPP: 領域を埋める
2016年08月29日(月)
なにがしかの計算に先立って device上の領域を特定の値で埋めておくって、よくやります。 0 で埋めるんだったら cudaMemset() がいちばんお手軽。
#include <cuda_runtime.h>
#include <iostream>

int main() {
  const unsigned int N = 8;
  int array[N];

  int* d_array;
  cudaMalloc(&d_array, N*sizeof(int));

  cudaMemset(d_array, 0, N*sizeof(int)); // 0で埋める

  cudaMemcpy(array, d_array, N*sizeof(int), cudaMemcpyDeviceToHost);

  for ( int item : array ) {
    std::cout << item << ' ';
  }
  std::cout << std::endl;

  cudaFree(d_array);
}
ですが困ったことにこの cudaMemset() はバイト単位で埋めるので、たとえばint配列のナカミを1で埋めることができません。上記コードの当該部分を
cudaMemset(d_array, 1, N*sizeof(int));
に書き換えると:



トホホな結果が得られます。

かといってお好みの値で埋めるのにこんなチンケなkernel書くのもアホくさい:
template<typename T>
__global__ void kernel_fill(T* to_fill, T value, unsigned int size) {
  unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
  if ( i < size ) to_fill[i] = value;
}
/*
   cudaMemset(d_array, 1, N*sizeof(int)); じゃダメだからかわりに:
   kernel_fill<<<(N+255)/256,256>>>(d_array, 1, N);
*/

NPP-libraryにおあつらえ向きのAPIが用意されてました:
#include <npps.h>
...
  int* d_array;
  cudaMalloc(&d_array, N*sizeof(int));

  // 32bit-signed int 領域: d_array[N] を 1 で埋める
  nppsSet_32s(1, d_array, N);


ここではint配列を埋めるべく nppsSet_32s() を使ってますが、
float なら nppsSet_32f, double なら nppsSet_64f などなど、
型に応じたAPIがひととおり揃ってます。

2016-08-29 10:40 | 記事へ | コメント(0) | トラックバック(1) |
| CUDA / NPP |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/126/

文字列の検索
2016年08月23日(火)
「長ぁい文字列:source の中から 文字列:target と一致する位置をぜんぶ見つける」ことを考えます。

CPU上でsingle-threadだとこんなカンジかしら:
  string source = "abracadabra";
  string target = "ra";
  vector<unsigned int> offsets; // 検索結果をココに求める

  for ( string::size_type pos = 0; 
        (pos = source.find(target, pos)) != string::npos;
        ++pos ) {
    offsets.push_back((unsigned int)pos);
  }
ここで使われている source.find(target, pos) は、source[pos]を検索開始位置としてtargetを探し、見つかった位置(なければstring::npos)を返します。

で、for-loopでは 見つかった位置+1 を新たな検索開始位置として次のtargetを探してます。

ですから n番目の検索は n-1番目の結果が得られていなくてはならんってことで、要するにマルチスレッド化しづらいコードなんですわ。

これをCUDAで求めるとなると、上記コードのような小賢しい(?)マネをせず、
愚直に「source+pos と target とが一致するか」を
pos = 0, 1 ... sourceの長さ について調べることにします。
そうすれば各threadが他のthreadの結果に依存しない、と。

__device__ unsigned int d_offsets_size = 0;

// src と dst が一致するか?
//    src[i] == dst[i] (i = 0, 1 ... size-1) なら true を返す
__device__ __forceinline__ bool equal(const unsigned char* src, 
                                      const unsigned char* dst, 
                                      unsigned int size) {
  for ( unsigned int i = 0; i < size; ++i ) {
    if ( *src != *dst ) { return false; }
    ++src;
    ++dst;
  }
  return true;
}

__global__ void kernel_findall(const void* source, unsigned int source_size,
                               const void* target, unsigned int target_size,
                               unsigned int* offsets) {
  unsigned int pos = blockDim.x * blockIdx.x + threadIdx.x;
  if ( pos <= source_size - target_size ) {
    // source+pos と target が一致するなら offsets に pos を追加する
    if ( equal((const unsigned char*)source + pos, (const unsigned char*)target, target_size) ) {
      offsets[atomicAdd(&d_offsets_size, 1)] = pos;
    }
  }
}

__host__ unsigned int launch_findall(const void* source, unsigned int source_size,
                                     const void* target, unsigned int target_size,
                                     unsigned int* offsets) {
  unsigned int result = 0;
  cudaMemcpyToSymbolAsync(d_offsets_size, &result, sizeof(unsigned int), 0);
  kernel_findall<<<(source_size+255)/256,256>>>(source, source_size, target, target_size, offsets);
  cudaMemcpyFromSymbol(&result, d_offsets_size, sizeof(unsigned int), 0);
  return result;
}

entry-modelグラボ:GT720(192 CUDA-core) でやってみました。
長さ10メガのsourceから 6文字のtarget を検索させたところ:



CPU single-thread で 20.7[ms] に対し GPUでは 12.2[ms] (内kernel本体は 8.6[ms]) となりました。 一昔前の安っすいグラボでもかなりがんばってくれてます。

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

cuRAND: device-APIでの乱数生成
2016年08月16日(火)
# 夏休みで間が空きましたゴメンナサイ

以前cuRANDによる乱数生成を紹介しました。
cuRAND host-API を使って大量の乱数を(device-codeなしで)吐かせました。

cuRANDの乱数生成にはdevice-APIてのもありまして、コチラはdevice-code上で各スレッド毎に乱数を生成します。

用意するのは乱数生成に必要なデータ:curandStateの列で、あらかじめcurand_init()で初期化しておきます。
__device__ ​ void curand_init(
                   unsigned long long seed, 
                   unsigned long long subsequence, 
                   unsigned long long offset,
                   curandState* state ) 

乱数の種をseed,初期状態をsubsequenceとoffsetで与えます。このふたつはスレッド毎に異なる値を与えておかないと、どのスレッドでも同じ乱数値が生成されます。
乱数生成アルゴリズムはXorshiftのバリエーション:XORWOWだそうな。

乱数生成は簡単:
__device__ ​ unsigned int curand ( curandState* state )
与えられたstateに基づいて32bit乱数を生成し、stateを更新します。他にも
- curand_uniform : 0〜1の一様分布
- curand_normal : 平均0, 標準偏差1 の正規分布
- curand_poisson : ポアソン分布
なんかが用意されてます。

おためしコード:十万個の一様乱数をdevie-APIで生成しました:

#include <device_launch_parameters.h>
#include <curand_kernel.h>

__global__ void kernel_init(curandState* state, unsigned long long seed, unsigned int size) { 
  unsigned int id = blockDim.x * blockIdx.x + threadIdx.x;
  if ( id < size ) {
    curand_init(seed, (unsigned long long)id, 0ULL, state+id);
  }
}

__global__ void kernel_generate_uniform(curandState* state, float* result, unsigned int size) {
  unsigned int id = blockDim.x * blockIdx.x + threadIdx.x;
  if ( id < size ) {
    result[id] = curand_uniform(state+id);
  }
}

// おためし: N個の一様乱数をdevice側で生成する
#include <cuda_runtime.h>
#include <iostream>
#include <vector>
#include <ctime> // time() for the seed
#include <cmath>

int main() {
  const unsigned int N = 100000;
  std::vector<float> random_seq(N);

  curandState* d_state;
  cudaMalloc(&d_state, N*sizeof(curandState));
  float* d_result;
  cudaMalloc(&d_result, N*sizeof(float));

  // d_stateを初期化し
  kernel_init<<<(N+255)/256,256>>>(d_state, time(nullptr), N);
  // 一様乱数を生成し
  kernel_generate_uniform<<<(N+255)/256,256>>>(d_state, d_result, N);
  // 結果(乱数列)をホストにコピーする
  cudaMemcpy(random_seq.data(), d_result, N*sizeof(float), cudaMemcpyDeviceToHost);

  // 最小値/最大値/平均値を求める
  float sum     =  0.0f;
  float min_val =  FLT_MAX;
  float max_val = -FLT_MAX;
  for ( float val : random_seq ) {
    sum += val;
    if ( min_val > val ) min_val = val;
    if ( max_val < val ) max_val = val;
  }
  std::cout << "N  : " << N       << std::endl 
            << "min: " << min_val << std::endl
            << "max: " << max_val << std::endl
            << "avg: " << sum/N   << std::endl;

  cudaFree(d_state);
  cudaFree(d_result);
  cudaDeviceReset();
}
Nsightで覗いたtimelineはこんなカンジ。



curand_init()にはそこそこの時間がかかりますが、curand_uniform()による乱数生成はすっごく速い、エントリ・モデルのGT720でも十万個で840マイクロ秒。device側でのシミュレーションなんかに向いてそうです。

2016-08-16 10:49 | 記事へ | コメント(0) | トラックバック(0) |
| CUDA / cuRAND |
トラックバックURL:http://blog.zaq.ne.jp/fareastprogramming/trackback/124/

NVTX: ユーザ定義のtimeline (おまけ)
2016年08月01日(月)
NVTX(NVIDIA Tools Extensions)によるオレオレtimelineはCUDAを使ってなくても有効ってわけでオアソビをひとつ。

NVTXの各関数はDLLで提供されています。.NETのDllImportアトリビュートを設定すればDLL内の関数をC#やVBから直接呼べるのを利用し、C#コードでtimelineを描かせてみよかと。
// NVTX.cs
using System.Runtime.InteropServices;

namespace NVTX {

  public static class Range {
    private const string nvtx_dll = "nvToolsExt64_1.dll";

    [DllImport(nvtx_dll, CallingConvention = CallingConvention.StdCall,
                         EntryPoint = "nvtxRangePushA")]
    private static extern int   nvtxRangePush([MarshalAs(UnmanagedType.LPStr)]string message);

    [DllImport(nvtx_dll, CallingConvention = CallingConvention.StdCall, 
                         EntryPoint = "nvtxRangePop")]
    private static extern int   nvtxRangePop();

    [DllImport(nvtx_dll, CallingConvention = CallingConvention.StdCall, 
                         EntryPoint = "nvtxRangeStartA")]
    private static extern ulong nvtxRangeStart([MarshalAs(UnmanagedType.LPStr)]string message);

    [DllImport(nvtx_dll, CallingConvention = CallingConvention.StdCall, 
                         EntryPoint = "nvtxRangeEnd")]
    private static extern void  nvtxRangeEnd(ulong id);

    [System.Diagnostics.ConditionalAttribute("NVTX")]
    public static void Push(string message) { nvtxRangePush(message); }
    [System.Diagnostics.ConditionalAttribute("NVTX")]
    public static void Pop() { nvtxRangePop(); }
    [System.Diagnostics.ConditionalAttribute("NVTX")]
    public static void Start(string message, ref ulong id) { id = nvtxRangeStart(message); }
    [System.Diagnostics.ConditionalAttribute("NVTX")]
    public static void End(ulong id) { nvtxRangeEnd(id); }
  }

}
こんな小さなinterfaceコードを用意し、フィボナッチ数列を求めるコードにtimelineを打ち込みます:
using NVTX;

namespace NVTX_marshal {

  class Program {
    static int Fib(int n) {
      NVTX.Range.Push(string.Format("Fib({0})",n));
      int result = ( n < 2 ) ? n : Fib(n-1) + Fib(n-2);
      NVTX.Range.Pop();
      return result;
    }

    static void Main() {
      ulong id = 0;
      NVTX.Range.Start("Main", ref id);
      System.Console.WriteLine("Fib(4) = {0}", Fib(5));
      NVTX.Range.End(id);
    }
  }

}
コンパイル・シンボル"NVTX"をセットしてコンパイルし、



nvvp上で起動すれば...



デキター♪

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

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/

次へ