krustf の雑記 Twitter

2012-02-22

CUDAでBoostを使ってみたかったんです

よくよく考えると、最近Boostをほとんどつかってないなーと思いました。ちゃんと言うと使いどころが出てこなかったんだと思います。

最近はCUDAをさわりつつ、という感じだったのでBoost.MPIとか使えると嬉しいなと思っています。なので少しずつBoost.MPIのドキュメント読んだり書いてみたりして使い方を覚えていく予定。

Host, Deviceの扱いが面倒ですね

CUDAってホストメモリとデバイスメモリと区別されていて、わりかしその辺気をつけないと簡単にエラーになったりとか、メモリ転送はDeviceToHost?HostToDevice?とか色々めんどくさいです。

見かけ上は、どちらもただのポインタとして振る舞うので非常に分かり辛い。そこはやはりC++、型を使っていい感じに扱いやすくしてほしいものです。

と言うわけで、結構適当な気もしますが書きました。Boostがうれしい。

環境は、nvcc 4.1, gcc 4.6.1 です。

#include <cstdio>
#include <boost/interprocess/smart_ptr/unique_ptr.hpp>
#include <boost/type_traits/arithmetic_traits.hpp>
#include <boost/type_traits/is_same.hpp>
#include <boost/static_assert.hpp>

namespace cuda { namespace detail {
  // detail used boost::interprocess.
  namespace ipc = boost::interprocess;

  // device, host delete.
  template<class T>
  struct host_deleter {
    void operator()(T* p) const {
      std::free(p);
    }
  };
  template<class T>
  struct device_deleter {
    void operator()(T* p) const {
      cudaFree(p);
    }
  };

  // device, host memory type getting.
  template<class T>
  struct memory_type {
    // T is arithmetic type.
    BOOST_STATIC_ASSERT(boost::is_arithmetic<T>::value);
    typedef ipc::unique_ptr<T, host_deleter<T> >   host_memory;
    typedef ipc::unique_ptr<T, device_deleter<T> > device_memory;
    typedef T element_type;
  };

  // unique_ptr index access.
  // allocate Host and Device memory.
  template<class T>
  typename memory_type<T>::host_memory make_hostmemory(std::size_t size) {
    // T is arithmetic type.
    BOOST_STATIC_ASSERT(boost::is_arithmetic<T>::value);
    return typename memory_type<T>::host_memory(
        reinterpret_cast<T*>(std::malloc(sizeof(T) * size))
    );
  }
  template<class T>
  typename memory_type<T>::device_memory make_devicememory(std::size_t size) {
    // T is arithmetic type.
    BOOST_STATIC_ASSERT(boost::is_arithmetic<T>::value);
    T* devp;
    cudaMalloc(reinterpret_cast<void**>(&devp), sizeof(T) * size);
    return typename memory_type<T>::device_memory(devp);
  }

  // memory data forward.
  template<class ElementType, class Dest, class Src>
  struct get_memcpy_kind {};
  template<class ElementType>
  struct get_memcpy_kind<ElementType
    , typename memory_type<ElementType>::host_memory
    , typename memory_type<ElementType>::device_memory
  > {
    static cudaMemcpyKind apply() { return cudaMemcpyDeviceToHost; };
  };
  template<class ElementType>
  struct get_memcpy_kind<ElementType
    , typename memory_type<ElementType>::device_memory
    , typename memory_type<ElementType>::host_memory
  > {
    static cudaMemcpyKind apply() { return cudaMemcpyHostToDevice; };
  };
  template<class ElementType>
  struct get_memcpy_kind<ElementType
    , typename memory_type<ElementType>::device_memory
    , typename memory_type<ElementType>::device_memory
  > {
    static cudaMemcpyKind apply() { return cudaMemcpyDeviceToDevice; };
  };
  template<class MemoryTypeL, class MemoryTypeR>
  void forward(MemoryTypeL const& dest, MemoryTypeR const& src, std::size_t size) {
    BOOST_STATIC_ASSERT((
      boost::is_same<typename MemoryTypeR::element_type,
                     typename MemoryTypeL::element_type>::value
    ));
    typedef typename MemoryTypeR::element_type element_type;
    typedef get_memcpy_kind<element_type, MemoryTypeL, MemoryTypeR> kind;
    cudaMemcpy(dest.get(), src.get(), sizeof(element_type) * size, kind::apply());
  }

  // return reference
  template<class HostMemoryType>
  typename HostMemoryType::element_type& at(HostMemoryType const& host, int index) {
    // T is host_memory.
    BOOST_STATIC_ASSERT((
      boost::is_same<HostMemoryType, 
        typename memory_type<typename HostMemoryType::element_type>::host_memory
      >::value)
    );
    return host.get()[index];
  }
  // return value
  template<class HostMemoryType>
  typename HostMemoryType::element_type at_c(HostMemoryType const& host, int index) {
    // T is host_memory.
    BOOST_STATIC_ASSERT((
      boost::is_same<HostMemoryType, 
        typename memory_type<typename HostMemoryType::element_type>::host_memory
      >::value)
    );
    return host.get()[index];
  }

}; // namespace detail

using detail::memory_type;
using detail::make_hostmemory;
using detail::make_devicememory;
using detail::at;
using detail::at_c;
using detail::forward;

}; // namespace cuda

template<class T>
__global__
void axpy_kernel(T a, T const* x, T const* y, T* z) {
  int tid = threadIdx.x;
  z[tid] = a * x[tid] + y[tid];
}

int main(int argc, char** argv) {
  static const int N = 11;

  // error! memory_type<T>, T is arithmetic type.
  // typedef cuda::memory_type<void>::host_memory void_host_memory;

  typedef cuda::memory_type<float>::host_memory   host_memory;
  typedef cuda::memory_type<float>::device_memory device_memory;

  host_memory x = cuda::make_hostmemory<float>(N);
  host_memory y = cuda::make_hostmemory<float>(N);
  for(int i = 0 ; i < N ; ++i) {
    cuda::at(x,i) = cuda::at(y,i) = i + 1;
  }

  device_memory devx = cuda::make_devicememory<float>(N);
  device_memory devy = cuda::make_devicememory<float>(N);
  cuda::forward(devx, x, N);
  cuda::forward(devy, y, N);

  // error! devx is device memory.
  // cuda::at_c(devx, 0);

  device_memory devz = cuda::make_devicememory<float>(N);
  axpy_kernel<<<1, N>>>(1.f, devx.get(), devy.get(), devz.get());

  host_memory z = cuda::make_hostmemory<float>(N);
  cuda::forward(z, devz, N);

  std::printf("result --- \n");
  for(int i = 0 ; i < N ; ++i) {
    std::printf("%f\n", cuda::at_c(z,i));
  }

  // deleted. z, devz, devy, devx, y, x
}

memory_typeというクラステンプレートでhost_memoryとdevice_memoryの型を取り出せます。typedef templateという奴ですか。

内部が若干「んー?」という感じになっているのはtypedefされたテンプレートの型を受け取って、そのtypedefされた型からテンプレートパラメータが欲しかったんですが、できないっぽかったので泣く泣くstatic assertなどで型チェックをしています。

C++11を使いたい所ですが、nvcc 4.1でもC++11の利用は難しいようなので、boost.interprocessからunique_ptrを引っ張ってホストメモリとデバイスメモリの型にしています。

この辺はうまいこと隠して、host_memoryとdevice_memoryはget()メンバ関数ぐらいしか使えない、コピーできない、とかドキュメントに書いておけば何の心配もないですね。

APIのベタ書きは一番嫌いな私からでした。

2012-02-18

Texture Memory使いたかったから勘弁してください

この記事書いてる途中にOpera for Ubuntu11.10が2度クラッシュしてデータ消えたのでもうUbuntuで記事書いてるときにOperaは使わないと心に決めました。

もう書くのめんどくさくなってきたのでさくっと書きます。

CUDAにはTexture MemoryがあってTexture Unitからのみアクセス可能ですが、Texture Unitには専用にキャッシュ領域があるので、アクセスが集中するようなデータをそこへ格納することでキャッシュの恩恵を受けてアクセスの高速化を図れるのでは、ということです。

例えば、行列ベクトル積 Ax = b においてベクトル x は計算時に複数回参照されることが分かります。そこにキャッシュの恩恵を適用すれば、アクセスにかかる時間が短くなり、計算が早くなるのではないか、というお話のようです。

データの正規化とか線形補間とかできるとかなんとかかかれていますが、私がお邪魔している研究室的には巨大な一次連立方程式を解くのが主目的なようなので、逆に補完されると困るって話なようなのでそういうのは使いません。

Texture Memoryを使用するには、textureクラステンプレートを宣言します。*1

  texture<float, 1, cudaReadModeElementType> float_tex;

パラメータが、テクスチャ内の数値型(int, int2, floatなど), 次元数, 読み取りモードの指定(cudaReadModeElementType, cudaReadModeNormalizedFloat)の順に指定します。

cudaReadModeNormalizedFloatが正規化を行うタイプで、私が使う予定はないです(多分...)

cudaReadModeElementTypeを指定しておくと、配列と似た感じのアクセスが可能になります。整数型のインデックスを指定したい場合、tex1Dfetch関数を使います。*2

次元数は、まあ1次元ベクトル2次元で行列、とか考えておくと分かるのではないでしょうか。

後はテクスチャにバインドしたりなんだったりなんですが、面倒なのでコード書けばいいですね。

static const std::size_t N = ...;

texture<float, cudaTextureType1D, cudaReadModeElementType> texture_;

template<class T>
__device__
void kernel_(int n, T const* x, T* y) {
  int i = ...;
  y[i] = x[i] * tex1Dfetch(texture_, i);
}

float *dev_vec; // デバイスのメモリ
float *dev_x, *dev_y;

cudaBindTexture(0, texture_, dev_vec,
                cudaCreateChannelDesc<float>(), sizeof(float) * N);
kernel_<<<..., ...>>>(N, dev_x, dev_y);
cudaUnbindTexture(texture_);

floatとかだったらこれで良いのですが、doubleはそのままでは渡せないようです。

CUDA 2.1 FAQ - NVIDIA Forums

このページのProgramming questions 23. "Can I read double precision floats from texture?"にTexture Memoryでのdoubleの使い方が書かれています。CUDA 2.1での解答なんですが、多分今もこの対応でいいのでしょう。

doubleを使おうとするとこうなります。

static const std::size_t N = ...;

texture<int2, cudaTextureType1D, cudaReadModeElementType> texture_;

template<class T>
__device__
void kernel_(int n, T const* x, T* y) {
  int i = ...;
  int2 v = tex1Dfetch(texture_, i);
  y[i] = x[i] * __hiloint2double(v.y, v.x);
}

double *dev_vec; // デバイスのメモリ
double *dev_x, *dev_y;

cudaBindTexture(0, texture_, dev_vec,
                cudaCreateChannelDesc<int2>(), sizeof(int2) * N);
kernel_<<<..., ...>>>(N, dev_x, dev_y);
cudaUnbindTexture(texture_);

アクセスめんどうですね。しかも型によって切り替えるし。研究室の人は多倍長整数なども実装したいと言っているらしく、多倍長にしたらまたアクセス方法が変わってきますね。

計算したい型情報から、テクスチャメモリの1要素のデータ型と、実際の計算で使うデータ型、後tex1Dfetchを使ってアクセスしてくれるようなTraitsクラスを書いた。

// デフォルト, テクスチャのデータ型と計算時のデータ型が同一
template<class DataType>
struct traits_ {
  typedef DataType value_type;
  typedef DataType texture_value_type;

  static const int fetch_dim = cudaTextureType1D;

  __device__
  value_type operator()(
    texture<texture_value_type, fetch_dim, cudaReadModeElementType>& vt,
    int index
  ) const {
    return tex1Dfetch(vt, index);
  }
};

// double 用
template<>
struct traits_<double> {
  typedef double value_type;
  typedef int2   texture_value_type;

  static const int fetch_dim = cudaTextureType1D;

  __device__
  value_type operator()(
    texture<texture_value_type, fetch_dim, cudaReadModeElementType>& vt,
    int index
  ) const {
    int2 v = tex1Dfetch(vt, index);
    return __hiloint2double(v.y, v.x);
  }
};

typedef double                    using_float_type;
typedef traits_<using_float_type> texture_traits;

texture<
    texture_traits::texture_float_type
  , texture_traits::fetch_dim
  , cudaReadModeElementType
> texture_;

static const std::size_t N = ...;

using_float_type *dev_vec;
using_float_type *dev_x, *dev_y;

cudaBindTexture( 0, texture_, dev_vec
               , cudaCreateChannelDesc<texture_traits::texture_value_type>()
               , sizeof(texture_traits::value_type) * N );
kernel_<<<..., ...>>>(N, dev_x, dev_y);
cudaUnbindTexture(texture_);

独自に定義した多倍長整数なんかでも、traitsクラスを特殊化すれば対応できると思う。

説明が結構端折ったり、足りてない部分があると思ったけど、CUDA C Programming Guide 4.1とか、ちゃんとドキュメント読まないと死ぬって事がよくわかってきた昨今です。

*1:CUDA C Programming Guideなどを読んでいると多分クラステンプレートという解釈で良いと思われます

*2:tex1DLayerdとかあるんだけどこれはなんだろうか

2012-02-17

Effective C++やC++11やテンプレートを教えてみた。

研究室院生の人に飲み会の席で冗談混じりに「教えて」とか言われたので、試しにやってみました。

CやJavaを学んだ+テンプレートは型Tをfloatやdoubleに置き換えて使っている程度ということだったので、まあEffective C++とかC++11とかテンプレートテクニック的な話だとかでいいかと思いつつそんなスライドを書きました。

ところどころ端折ってたり、説明が足りない、間違っているがあるかも知れない。多分大丈夫。色付けはいくつか忘れている。

[誤字報告] 2012/02/18

  1. "default and deleted function" → "defaulted and deleted functions"

ideone C++0xモード便利です。本当に。

今回説明しなかったけど説明しといたほうが良かった箇所

  1. explicit = スライドには説明ないけど口頭+ホワイトボードで行った。
  2. constexpr
  3. move semantics, rvalue reference

2012-02-03

メモ:OpenGL Interoperability の話

長いので Interop にして話す。OpenGLとの相互運用のこと。Direct3Dとも可能だが、Direct3Dはいかんせん初期化が面倒なので使いたくないので無視。

CUDA C Programming Guideの3.2.11辺りにGraphics Interoperabilityのページがある*1のでそれと、SDKのサンプルを見ると分かりやすい。後CUDA by Exampleとか。大学図書館に英語の原書があったので、借りて英語がよく分からないながらも読んでます。

正直、SDKのサンプルってDirectXとかもそうだけどなんかコード自体に読ませる気がない長さを持っている気がする。どうしてあんなに長いのか。サンプル名に simple などかかれていても全くシンプルに見えないので困る。

まあ、ちょっと長いほうがコードを読む甲斐があるような気もしますが、そう思わないとやってられない。

で、OpenGL初期化とか描画周り差し引けばそこまで難しくない。むしろ描画周りのパラメータの設定とかのほうがしんどいと思う。

まず、CUDAの関数OpenGL初期化よりも初めに、cudaGLSetGLDeviceで使用するデバイスを決定する。

cudaGLSetGLDevice(0);

CUDA by Exampleを読むと、compute capabilityを指定して対応するデバイスのIDを取得したりしていたが、まあとりあえずそこまでする必要はないと思うので省略。0 で多分デフォルトデバイスになるはず。

その後、VBO(Vertex Buffer Object)を作って、それをカーネルで使えるように登録する。

int width  = ...; // VBO のサイズ
int height = ...;

GLuint vbo;
struct cudaGraphicsResource* vbo_cuda;

glGenBuffer(1, &vbo);
glBindBuffer(GL_ARRAY_BUFFER, &vbo);

unsigned int size = width * height * 4 * sizeof(float);
glBufferData(GL_ARRAY_BUFFER, size, 0, GL_DYNAMIC_DRAW);
glBindBuffer(GL_ARRAY_BUFFER, 0);

// カーネルで使えるように登録
cudaGraphicsGLRegisterBuffer( &vbo_cuda, vbo
                      , cudaGraphicsMapFlagsWriteDiscard );

で、カーネルへ渡すときはリソースのビューを作成する。

float4* v;
cudaGraphicsMapResources(1, &vbo_cuda, 0);
size_t bytes;
cudaGraphicsResourceGetMappedPointer((void**)&v, &bytes, vbo_cuda);

この v がVBOへの先頭のポインタとなるので、これに対してカーネルで頂点情報を書き込めば良い。

書き込んだ後はUnmapしないと描画にも使えないので必ず行う。

kernel<<< ... , ... >>>(v, width, height, ...);

cudaGraphicsUnmapResources(1, &vbo_cuda, 0);

後は、VBOを使って正しく描画するだけ。アプリケーション終了時に必ずリソースを破棄するようにstd::atexitを使って終了関数を登録しておく。

void cleanup() {
    // 登録の解除
    cudaGraphicsUnregisterResource(vbo_cuda);
    glBindBuffer(GL_ARRAY_BUFFER, vbo);
    glDeleteBuffers(1, &vbo);
}

std::atexit(cleanup);

という感じでいいらしい。Interop よりも綺麗に描画しようとするとシェーダとか使わないといけなくなるんでしょうねえ。simpleGLしか読んでなかったけど、oceanFFTとか読んでみますか。

*1:version 4.0での話

2012-02-02

メモ:neocomplcacheでC++のインクルード補完

なぜかすごい躓いたのでとりあえずメモ。

.vimrcでg:neocomplcache_include_paths.cppに設定する方法もあるようだが、何かうまくいかないのでとりあえずpathに追加する形で。

" ~/local/includeにboostのインクルードディレクトリが置いてある
set path+=$HOME/local/include
" g++が参照している標準ライブラリのヘッダパス
" 'gcc -v' で --with-gxx-include-dir のそれ
set path+=/usr/include/c++/4.6

高校から使わされて使い始めてからvimプラグインなんて使わなくても生きていけてた人だったので全然慣れてませんが、少しずつ慣れていきましょう。