2016年8月8日月曜日

GpuMatの内部を探検してみる

はじめに

OpenCVにはGpuMatというCUDA実装を行うためのデータ構造が用意されており、CUDAを使って実装された各種アルゴリズムもcudaモジュールという形で提供されています(※詳細は公式ドキュメントを参照ください)。

この記事ではGpuMat内部で行われている処理とGpuMatを使う上でのTipsをいくつか紹介します。特に

  • OpenCVにあるCUDA実装を読めるようになりたい
  • GpuMatと自作CUDAコードを組み合わせたい

という方には参考になるかもしれません。

前置き

まずは、cudaimgprocモジュールの

cv::cuda::cvtColor(d_src, d_dst, cv::COLOR_BGR2GRAY);

を呼ぶだけの単純なプログラムを実行して、Nsightで見てみましょう。
すると、以下のようなCUDAカーネルが動いていることがわかります。



これを見るとtransformSimpleGlobPtrという見覚えのない文字列が出てきて戸惑うかもしれませんね(私は初見でかなり戸惑いました・・・)。ただ、この記事を読み終わる頃にはこれらのパラメータの意味がわかるようになっているはずです。

インスタンス生成

GpuMatクラスのインスタンス生成でよく使う方法を紹介します。

以下はrows、cols、typeの順に指定する例です(この初期化方法では、width、heightの順ではなく、rows、colsの順である点に注意が必要です)。

cv::cuda::GpuMat d_img(240, 320, CV_8UC1);

次にサイズ(width、height)、typeの順に指定する例です。

cv::cuda::GpuMat d_img(cv::Size(320, 240), CV_8UC1);

次にホスト側のimgと同じ画像バッファの値を持ったGpuMatのインスタンスを生成する例です。

cv::Mat img(cv::Size(320, 240), CV_8UC1);
cv::cuda::GpuMat d_img(img);

データ転送

ホスト、デバイス間のデータ転送はGpuMatクラスのupload、downloadメソッドを用います。uploadメソッドがホスト→デバイスの転送、uploadメソッドがデバイス→ホストの転送を行うメソッドとなっています。

cv::Mat img(cv::Size(320, 240), CV_8UC1);
cv::cuda::GpuMat d_img(src1.size(), src1.type());

d_img.upload(img);   //ホストからデバイスに転送
d_img.download(img); //デバイスからホストに転送

download、uploadメソッドは、内部的にはcudaMemcpy2Dを呼んでいるだけです。ただし、download、uploadメソッドの引数としてstreamを指定した場合はcudaMemcpy2DAsyncが呼ばれます。

cudaモジュールの内部実装

全ての実装がこうなっているというわけではありませんが、cudaモジュールで提供されるAPIの内部処理をおおまかにまとめると以下のような感じになっています。



プリミティブな処理に関しては何気にNPP(NVIDIA Performance Primitives library)に投げていたりします。また、この図に突然出てきたgridTransformについては後ほど述べることとします。そのため、この時点ではそういうものがあるんだなという理解でOKです。


cudevモジュール

CUDAカーネルを書く上でよく用いるデータ構造、処理などはcudevモジュールにまとまっています。
ここではその中でも重要なものについてピックアップして紹介します。

GlobPtrSz構造体

CUDAカーネルの引数として画像バッファ、width、height、strideなどを渡すことが多いと思います。
これらを一つ一つカーネル引数で渡すとカーネル引数が煩雑になりがちです。

cudevモジュールではこれらのパラメータをまとめたGlobPtrSz構造体が提供しています。
文章だけだとピンとこないと思いますのでサンプルコードを交えながら紹介します。

GlobPtrSz構造体を生成するサンプルコード

cv::cudev::globPtrメソッドに画像バッファのポインタ、step、rows、colsを渡すと、cv::cudev::GlobPtrSz構造体のデータを作ることができます。GlobPtrSz構造体を生成し、CUDAカーネルに渡すサンプルコードを以下に示します。

void launchMyKernel(cv::cuda::GpuMat &src, cv::cuda::GpuMat &dst)
{
    cv::cudev::GlobPtrSz<uchar> pSrc = 
        cv::cudev::globPtr(src.ptr<uchar>(0), src.step, src.rows, src.cols * src.channels());
    cv::cudev::GlobPtrSz<uchar> pDst = 
        cv::cudev::globPtr(dst.ptr<uchar>(0), dst.step, dst.rows, dst.cols * dst.channels());

    const dim3 block(64, 2);
    const dim3 grid(cv::cudev::divUp(src.cols, block.x), cv::cudev::divUp(src.rows, block.y));

    myKernel<<<grid, block>>>(pSrc, pDst);

    CV_CUDEV_SAFE_CALL(cudaGetLastError());
    CV_CUDEV_SAFE_CALL(cudaDeviceSynchronize());
}

上記サンプルコードのようにCUDAカーネル(myKernel)の引数としてcv::cudev::GlobPtrSz構造体の変数を渡すことでCUDAカーネルの引数がシンプルにできていることがわかります。

GlobPtrSz構造体を参照するサンプルコード

cv::cudev::GlobPtrSz構造体のメンバ変数data、step、rows、colsを参照することで必要なデータにアクセスすることができます。GlobPtrSz構造体を参照するサンプルコードを以下に示します。

__global__ void myKernel(const cv::cudev::GlobPtrSz<uchar> src, cv::cudev::GlobPtrSz<uchar> dst)
{
    int x = blockDim.x * blockIdx.x + threadIdx.x;
    int y = blockDim.y * blockIdx.y + threadIdx.y;
    dst.data[y*src.step + x] =  UCHAR_MAX - src.data[y*src.step + x];
}

GpuMatにある画像バッファのアドレス取得

GpuMatクラスのインスタンスは画像バッファ(デバイス側)のアドレスを保持しています。GpuMatクラスにある画像バッファのアドレス取得方法について代表的な方法は以下の3パターンあります(これ以外で別の方法をご存知の方がいらっしゃったらぜひ教えてください)。
  1. メンバ変数dataから取得
  2. ptrメソッドから取得
  3. GlobPtrSz構造体のメンバ変数dataから取得
以下にこれらの取得方法を行ったサンプルコードを示します。
cv::cuda::GpuMat src(cv::Size(320, 240), CV_8UC1);

// (1) メンバ変数dataから取得
uchar *pData1 = src.data;

// (2) ptrメソッドから取得
uchar *pData2 = src.ptr<uchar>();

// (3) GlobPtrSz構造体のメンバ変数dataから取得
cv::cudev::GlobPtrSz<uchar> pSrc = 
 cv::cudev::globPtr(src.ptr<uchar>(), src.step, src.rows, src.cols * src.channels());
uchar *pData3 = pSrc.data;

gridTransform

ここでは、これまでの説明で先延ばしにしてきたgridTransformについて紹介します。 cudevモジュールにて定義されるgridTransformは大別すると以下の3つに分類されます。
  • gridTransformUnary
  • src、dst、opを引数に持ち、「入力1つ、出力は1つ。入力に対してopを適用」といったケースのカーネルを呼ぶのに用いられる
  • gridTransformBinary
  • src1、src2、dst、opを引数に持ち、「入力2つ、出力は1つ。入力に対してopを適用」といったケースのカーネルを呼ぶのに用いられる
  • gridTransformTuple
この記事の冒頭でNsightの結果に出てきたtransformSimpleは、src、dst、op、maskといったパラメータを持ち、様々な処理を統一的に呼べるよう汎用性を持たせたCUDAカーネルとなっています(そのため、Nsightの結果にtransformSimpleが表示されていたわけですね)。 cv::cuda::cvtColor(BGR2GRAY)を例に挙げると、厳密には以下のようなシーケンスで処理が行われます。
cv::cuda::cvtColor
 BGR_to_GRAY
  cv::cuda::device::BGR_to_GRAY_8u
   cv::cudev::gridTransformUnary
    cv::cudev::grid_transform_detail::transform_unary
     cv::cudev::grid_transform_detail::TransformDispatcher
      cv::cudev::grid_transform_detail::TransformSimple
       →ここでCUDAカーネルが呼ばれる
grid_transform_detail以下の振る舞いは長くなりそうなのでここでは詳細について割愛します。 また、cv::cudev::grid_transform_detail::TransformSimpleの引数opに目的の処理(上記例だとBGR2GRAYの色変換処理)の関数ポインタが渡されているため、transformSimpleでは目的の処理に対応したCUDAカーネルが呼ばれるようになっています。 このBGR2GRAYの色変換の例ではmodules/cudev/include/opencv2/cudev/functional/detail/color_cvt.hppにて定義される以下の処理が呼ばれるようになっています。
template <typename T, int scn, int bidx> struct RGB2Gray
        : unary_function<typename MakeVec<T, scn>::type, T>
{
    __device__ T operator ()(const typename MakeVec<T, scn>::type& src) const
    {
        const int b = bidx == 0 ? src.x : src.z;
        const int g = src.y;
        const int r = bidx == 0 ? src.z : src.x;
        return (T) CV_CUDEV_DESCALE(b * B2Y + g * G2Y + r * R2Y, yuv_shift);
    }
};

アトミック関数

CUDAが提供するatomic関数はCompute Capabilityによってfloat型をサポートしていなかったりするのでatomicCASを使ってatomic関数を自作したことがある方もいらっしゃるかもしれません。OpenCVのcudevモジュールではCompute Capabilityの違いを吸収した以下のアトミック関数を提供しています。
  • cv::cudev::atomicAdd
  • cv::cudev::atomicMin
  • cv::cudev::atomicMax
そのため、上記APIを使うことでユーザがatomicCASを使ってアトミック関数自作する必要がなくなります。とはいえ、Compute Capabilityが2.0以降であれば普通にfloat型サポートしているので、今となってはこれらの関数の意義は大分薄くなってる気もします・・・。

Tips

ここでは個人的にまとめていたGpuMatを使う上でのTipsを紹介します。

GpuMatを入力にした自作CUDAカーネルを作りたい

GpuMatを入力として自作CUDAカーネルを動かす単純なサンプルコードを下記リポジトリに置いていますのでご参考ください。 https://github.com/atinfinity/launchMyKernel 自作CUDAカーネルはネガポジ反転するだけの非常にシンプルなものですが、この方法を応用すると「ある処理はOpenCV(cudaモジュール)に任せて、ある処理は自前で書いてチューニングする」みたいなことができるようになります。

OpenCVのCUDAカーネル内をデバッグしたい

CMakeでWITH_CUDA=ONとしてOpenCVをビルドした場合、通常、NVCCでOpenCVのCUDAカーネルをビルドする際にデバッグ情報が付与されないのでOpenCVのCUDAカーネル内をデバッグできません。 デバッグできるようにする最も簡単な方法はcmake/OpenCVDetectCUDA.cmakeを編集し、NVCCのコンパイラオプションを追加するという方法です。具体的には、
# NVCC flags to be set
set(NVCC_FLAGS_EXTRA "")
となっている箇所を
# NVCC flags to be set
set(NVCC_FLAGS_EXTRA "-G -g")
に変更した後、CMakeを実行し、OpenCVをビルドすることでOpenCVのCUDAカーネル内をデバッグできるようになります。

複数GPUを使う

cv::cuda::setDevice関数を使うことでデバイス指定することができます。 また、公式サンプルにも複数GPUを使うものがあるので参考になるかもしれません。
これらのサンプルは、TBBやpthread等を用いて処理をデバイス毎に振り分けるシンプルなものになっています。

デバイスのfeatureサポートをチェックする

動作デバイスのfeatureサポートをチェックするためのcv::cuda::deviceSupports関数が提供されています。チェックできるfeatureは公式ドキュメントにまとまっています。以下にWarpShuffleが使えるかどうかをチェックするサンプルコードを紹介します。
bool isSupportWarpShuffle = 
    cv::cuda::deviceSupports(cv::cuda::FeatureSet::WARP_SHUFFLE_FUNCTIONS);

GpuMatの画像データのウィンドウ表示でちょっと楽をする

GpuMatのデータをhighguiのウィンドウに表示する際、通常、downloadメソッドでホストに転送して、cv::imshowで表示すると思いますがちょっとだけ簡単にできる方法があります。 具体的には以下のような設定を行うことでユーザが明示的にホストへ転送しなくてもよくなります。
  • OpenCVをビルドする際にWITH_OPENGLを有効にする
  • cv::namedWindow関数にセットするflagにcv::WINDOW_OPENGLを付与する
以下にサンプルコードを示します。
#include <opencv2/core.hpp>
#include <opencv2/core/cuda.hpp>
#include <opencv2/cudaimgproc.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/highgui.hpp>

#include <iostream>

int main(int argc, char *argv[])
{
    cv::Mat src = cv::imread("lena.jpg", cv::IMREAD_UNCHANGED);
    cv::Mat dst;

    // 画像の読み込みに失敗したらエラー終了する
    if(src.empty())
    {
        return -1;
    }

    cv::cuda::GpuMat d_src(src), d_dst;
    cv::cuda::cvtColor(d_src, d_dst, cv::COLOR_BGR2GRAY);

    // highgui(normal)を使った表示
    d_dst.download(dst); // ホストメモリに転送する
    cv::namedWindow("normal", cv::WINDOW_AUTOSIZE);
    cv::imshow("normal", dst);

    // highgui(OpenGL)を使った表示
    // この方法だと明示的にホストに転送する必要がない
    cv::namedWindow("highgui(OpenGL)", cv::WINDOW_AUTOSIZE | cv::WINDOW_OPENGL);
    cv::imshow("highgui(OpenGL)", d_dst);

    cv::waitKey(0);
    cv::destroyAllWindows();

    return 0;
}
詳細については、筆者の別記事にて解説していますのでご参考ください。

cudaモジュールのビルド時間を短縮する

OpenCVのcudaモジュールをビルドしたことがある方は経験があるかもしれませんが、普通にCMakeでWITH_CUDAをONにしてビルドすると非常にビルド時間が掛かってしまいます(マシンスペックによりますが、以前6時間以上掛かったことも・・・)。これはOpenCVビルド時に複数のCompute Capability向けにビルドをしているというのが理由の一つです。 そのため、CMakeオプションにて対象デバイスのCapability番号を明示的に指定する(例えばCUDA_ARCH_BIN="3.5"と指定する)ことで、特定のCompute Capability向けのコンパイルにとどめることができ、ビルド時間を大きく短縮することができます。お手持ちのNVIDIA GPUに対応するCompute CapabilityはNVIDIAサイト調べることができます。 また、必要なcudaモジュールだけビルドすることでもビルド時間を短縮することができます (例えば、CMakeオプションBUILD_opencv_cudaarithmはONにして、BUILD_opencv_cudaoptflowはOFFにするなど)。 こちらでOpenCVのビルドスクリプトを公開しているのでご参考ください。

おわりに

この記事ではGpuMatの内部でどんな処理が走っているかとGpuMatを使う上でのTipsをいくつか紹介しました。この記事が、OpenCVにあるCUDA実装を読み解いたり、GpuMatと自作CUDAコードを組み合わせる際の助けになりましたら幸いです。

2016年8月2日火曜日

CUDAをRadeonで動かす(導入編)

みなさん、今日も元気にGPGPUしていますか?

去年(SC15)の話ですが、「RadeonでCUDAが使えるようにするよ!」とAMDが発表したニュースを覚えている方いらっしゃいますでしょうか。Boltzmann Initiativeの話です。

あれからしばらく時が流れ、機は熟しました。「CUDAはNVIDIA専用」そんな時期は今は昔。そう、CUDAをAMDのGPUであるRadeonで動かすことに成功しました!!

ので、ここではその方法と、その時に使ったコードを紹介したいと思います。

予めお断りしておくと、表題にもある通り、今回は導入編ということで「動くことを確認する」までになります。 そのうち時間を見つけて「性能測定編」もやってみたいと思いますが、ぜひこれを読んだ方、ぜひ性能測定してみてください&結果を教えて下さい!

また、この情報は、本記事執筆時点(つまり2016年7月21日現在)の話です。 Boltzmann Initiativeは現在も活発に進行中で、未来では手順が変わったりする可能性が十分にありますので、ご了承ください。

参考文献

説明に入る前に、最初に参考文献を挙げておきます。 本記事はこれらの文献のまとめになりますので、手順書さえあればできるという方・最新かつ正確な情報がほしい方は、以下を参照してください。

概要

使うのはROCmというプロジェクトです。 ROCmとは、Radeon Open Computeの略のようです。

ROCmはいくつかのプロジェクトのプロジェクトの集合体となっていますが、今回直接的に使うのは以下の3つです

  • Kernel Fusion DriverおよびRuntime群
  • HCC compiler
  • HIP

Kernel Fusion Driver

Kernel Fusion Driver(KFD)とは、ROCmの各種プロジェクトを動かすための「Linuxカーネルと一体となったドライバー」のことです。

なので、後述の手順で説明しますが、カーネルを入れ替える必要があります。 個人的にはkmodでどうにかしてほしいんですが・・・。

ともかく、このKFDを入れることで、AMD Radeonの上で様々なことができるようになります。

HCC compiler

HCC compilerは、KFDの上で動く、異種混合実行環境向けのC++コンパイラーです。 HCCとは、Heterogeneous Compute C++の略のようです。

現在のところ、バックエンドとしてはAMD独自のAMD GCN ISAを生で出力するモードと、HSA財団の策定したHSAILという中間言語を出力するモードがあります。 ただし、今回はRadeonを使うので前者で十分なことと、既定で前者が使われること、加えてそもそも後者はこの先開発はされないようなので、ここではAMD GCN ISAを使うモードを前提として話を進めます (とは言っても、バックエンドが何であるかはアプリケーション層ではあまり気にする必要はなさそうですので、よく分からない方は飛ばしても大丈夫です)

また、フロントエンドにはC++17 Parallel STL、C++ AMP、OpenMP offload、HC C++ APIなどが用意されています。 今回は説明しませんが、C++ AMPをLinuxで動かしたい場合は、このHCC compilerを使うと動かせるようになります。

HIP

HIPは、CUDAをHCCでコンパイルできるようにするプラットフォームです。今回の話の主役になります。何の略称かは分かりません・・・。

HIP自体は、ほぼCUDAと同じ形式のRuntime APIやカーネル言語を既定するもので、基本的にcuXXXXをhipXXXXに置き換えたものです(例えばcuMemcpyhipMemcpyに対応します)。

HIPで書かれたものは、hipccというコンパイラによって解釈され、バックエンドであるnvccまたはhccに渡すことで、NVIDIAでもAMDでも動かせるようにできます。 このhipccは、CUDAあるいはHC言語そのものへの直接的な翻訳機に近い動きをするので、HIPで書くことにオーバーヘッドなどはなく、性能劣化することは全くないらしいです。

また、hipifyという、CUDAからHIPへ(逆)変換するツールが付属しており、これを用いることで既存のCUDAコードを簡単にHIPが解釈できるようになります。

ということで、このHIPを信じるなら、CUDAで書いたコードを、性能劣化させることなくRadeonでも動く形にすることができる、ということになります。

手順

ということで、実際に動かしてみましょう。

環境の準備

ROCmのREADMEROCm Hardware Requirementsに情報が散乱しているのですが、ROCmを動かすには、現在のところ以下の環境が必須です。

OS
Ubuntu 14.04.04またはFedora23
(おそらく動く:Ubuntu 16.04またはFedora 22)
CPU
Intel Haswell以降またはAMD Zen以降(PCIe Gen3 Atomicsのサポート必須)
GPU
Fiji以降のAMD GPU

今回は、Ubuntu 14.04、Intel Core i7-4790、AMD Radeon R9 FuryXの環境で試しました。

また、カーネルを入れ替えるという性質上、クリーンインストールすることにしました。 既存のシステムの上に入れても動かないわけではないと思いますが、何かあった時にサルベージが大変なので・・・。

あとこれ重要なのですが既存のAMDのプロプライエタリ・ドライバーとは衝突してまともに動かなくなるっぽいです。 公式ブログのコメント欄でしれっと書かれているのですが、 KFDでドライバーも提供されているので、ROCm導入時には既存のプロプライエタリ・ドライバーは除去する必要があります。 つまり、ROCmを入れるとOpenCLは使えなくなるということです(Issue参照)。 将来的にはサポートされるかもしれませんが・・・。

ということで、既存の環境をかなり破壊することになるので、試される方は新しいストレージを持ってきて、クリーンインストールして始めることをオススメします (これで安全にはなりますが、カーネルの入れ替え等に限らず、ここに書いてあることは全て自己責任でお願いします)

ROCmのインストール

計算機を用意して、OSをクリーンインストールしたら、まずは普通に起動します。

起動したら、とりあえずまずインストールされているパッケージを最新に更新します。 これはセキュリティーリスクもそうですが、今回用意した環境だと後でバージョン違いでやり直しという羽目になったので、ここで一度更新しておくことをオススメします。 Ubuntuならapt-get upgrade、Fedoraならdnf updateですね。

その後、ビルド済みバイナリをAMDのレポジトリから持ってきます。 Ubuntuの場合Fedoraの場合の両方が用意されています。 今回はUbuntuを使ったので、以下のコマンドを実行しました。

wget -qO - http://packages.amd.com/rocm/apt/debian/rocm.gpg.key | sudo apt-key add -
sudo sh -c 'echo deb [arch=amd64] http://packages.amd.com/rocm/apt/debian/ trusty main > /etc/apt/sources.list.d/rocm.list'
sudo apt-get update
sudo apt-get install rocm

なお、AMDのサーバーがとても遅いのか、恐ろしく時間がかかります。 具体的には、一晩放置しても終わってませんでした。半日まるっとかかるぐらいは覚悟して放置しましょう。

インストールが無事に終わったら再起動します。 再起動後、uname -rとかでカーネル名が4.4.0-kfd-compute-rocm-rel-1.1.1-10などになっていれば新しいカーネルに切り替わっています。

その後、手順に従って、HSAILが動くことを確認しましょう。 ちゃんと動けばインストール成功です。 これで、ROCmの一式が全て入っているので、HCC compilerもHIPも全て/opt/rocmの下にインストールされました。

hipify-clangの導入

先のレポジトリから持ってくる手順で、CUDAからHIPに変換するツールであるhipifyは/opt/rocm/bin/hipifyにインストールされました。 が、実はこの標準のhipifyだけでは変換が不十分になります。 これは、このhipifyがただの文字列置換Perlスクリプトであり、C++の構文解析を一切していないからです。

なので、hipify-clangと言う、clangベースのツールを別に入れる必要があります。 別のツールといっても、clang-hipifyもHIPの中の1つです(↑の手順で入らないのはまだβ版だからのようですね)。

手順は、HIPのmasterレポジトリをダウンロードして、 clang-hipifyのREADMEに従ってビルド&インストールしましょう。 今回は以下のコマンドを実行しました。

cd ~
git clone https://github.com/GPUOpen-ProfessionalCompute-Tools/HIP.git

cd HIP
wget http://llvm.org/releases/3.8.0/clang+llvm-3.8.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz
tar xvfJ clang+llvm-3.8.0-x86_64-linux-gnu-ubuntu-14.04.tar.xz

mkdir build
cd build
cmake -DBUILD_CLANG_HIPIFY=1 -DLLVM_DIR=~/HIP/clang+llvm-3.8.0-x86_64-linux-gnu-ubuntu-14.04/ -DCMAKE_BUILD_TYPE=Release ..
make
sudo make install

wget http://developer.download.nvidia.com/compute/cuda/repos/ubuntu1404/x86_64/cuda-repo-ubuntu1404_7.5-18_amd64.deb
sudo gdebi cuda-repo-ubuntu1404_7.5-18_amd64.deb
sudo apt-get update
sudo apt-get install cuda-minimal-build-7-5 cuda-curand-dev-7-5

これで、/opt/rocm/hip/bin/clang-hipifyがインストールされます。

元にするCUDAコードの用意

ということで準備が終わったので、実際にCUDAコードを動かす手順に移ります。

ということで手始めに元にするコードを用意します。今回は簡単なベクトルの加算を用意しました。

__global__
void kernel(double z[], const double x[], const double y[])
{
 const auto i = (blockIdx.x * blockDim.x + threadIdx.x);

 z[i] = x[i] + y[i];
}

せっかくなので、これをまずは普通にCUDAの動く環境で動かしてみてください(vectorAddに移動してmakeで動くはずです)。

CUDAをHIPに変換する

CUDAコードを用意したら、いよいよ次はCUDAをHIPに変換する時です!

ということで、とりあえず普通についてくるhipifyを使ってみましょう。 hipify [yoursource.code]とやると、標準出力に変換結果が表示されます。 これをやると、例えば

#include "kernel.hpp"

__global__
void kernel(double z[], const double x[], const double y[])
{
 const auto i = (blockIdx.x * blockDim.x + threadIdx.x);

 z[i] = x[i] + y[i];
}


void Kernel(double z[], const double x[], const double y[],
        const std::size_t block, const std::size_t thread)
{
 kernel<<<block, thread>>>(z, x, y);
}

は、

#include <hip/hip_runtime.h>
#include "kernel.hpp"

__global__
void kernel(double z[], const double x[], const double y[])
{
        const auto i = (hipBlockIdx_x * hipBlockDim_x + hipThreadIdx_x);

        z[i] = x[i] + y[i];
}


void Kernel(double z[], const double x[], const double y[],
        const std::size_t block, const std::size_t thread)
{
        hipLaunchKernel(HIP_KERNEL_NAME(kernel), dim3(block), dim3(thread), 0, 0, z, x, y);
}

と変換されます。概ね一対一になっていて対応関係はそんなに難しいものではないと思います。

ただ、これは実は不完全です。これを後述のhipccに食わせてもコンパイルできません。 やってみると分かるのですが、カーネル関数の第一引数に必要な物がない、と言われます。 実はHIPの言語仕様として、カーネル関数の最初にはhipLaunchParm lpというものが必要ということになっています。

そのため、通常のhipifyを使っただけでは、変換後に手動でhipLaunchParm lpを挿入する必要があります。面倒ですね。

そこで、先にインストールしたhipify-clangの出番です。 これを使ってhipify-clang [yoursource.code] -o [yourhip.code]とすると、無事に全部変換してくれます。

しかしhipify-clangだけだとまだ不完全なようで、ひとつは、hip_runtime.hのincludeを追加してくれないんですね・・・。 ということで、現時点では正しく自動変換させるには、hipify-clangを書けた後に、通常のhipifyをする必要があるようです。

HIPのビルド

HIPに変換できたら、あとはhipcc [yourhip.code]とするだけでビルドできます。 あとはビルドされた実行ファイルを普通に実行できるはずです。

というわけで、HIPへの変換からビルドまでを一括でやってくれるMakefileを書きました。 ご覧のとおり、このMakefileでは以下のことをやってくれます。

  1. カレントディレクトリの*.cppをHIP Runtime APIを使うように変更し、obj/*.hip.cppに保存
  2. カレントディレクトリの*.cuをHIP Kernel Languageを使うように変更し、obj/*.hip.cuに保存
  3. 各obj/*.hip.c(pp|u)を、hipccでビルドしてobj/*.c(pp|u).oにコンパイル
  4. 全てのobj/*.c(pp|u).oをリンクして実行ファイルを出力

先述のベクトル加算の例にこのMakefileを置いたものも用意しました。 makeとやると、vectorAddという実行ファイルができます。

$ git clone git@bitbucket.org:LWisteria/hipvectoradd.git
$ cd hipvectoradd/vectorAdd

$ make
/opt/rocm/hip/bin/hipify-clang main.cpp -o main.cpp.hip
hipify main.cpp.hip -inplace
rm main.cpp.hip.prehip -f
mv main.cpp.hip obj/main.hip.cpp
hipcc -g -O3 -Wall -I/opt/rocm/include/hip -std=c++11 -I. -c obj/main.hip.cpp -o obj/main.cpp.o
/opt/rocm/hip/bin/hipify-clang kernel.cu -o kernel.cu.hip
hipify kernel.cu.hip -inplace
rm kernel.cu.hip.prehip -f
mv kernel.cu.hip obj/kernel.hip.cu
hipcc -g -O3 -Wall -I/opt/rocm/include/hip -std=c++11 -I. -c obj/kernel.hip.cu -o obj/kernel.cu.o
hipcc -o vectorAdd obj/main.cpp.o obj/kernel.cu.o

$ ./vectorAdd
Device: Fiji
0: -1.38309, -1.38309
1: -1.91112, -1.91112
2: -0.682073, -0.682073
3: 1.17407, 1.17407
4: 1.95078, 1.95078
5: 0.93395, 0.93395
6: -0.941547, -0.941547
7: -1.95139, -1.95139
8: -1.16713, -1.16713
9: 0.69018, 0.69018
  

実行結果はあまりちゃんと書いてませんが、main.cppを読んでもらうと分かる通り、左側の数値がCPUの実行結果で、右側がGPUでの実行結果です。 ちゃんと実行できたことが分かりますね!

おわりに

ということで、CUDAのコードをRadeonで動かすことに成功しました!

ROCmあるいはHIPに関しての日本語での情報がほぼ皆無だったので、この記事が誰かの役に立てばとても嬉しく思います。

なお、最初に書いた通り、今回は「動かす」ことを目的にしたので、特に性能評価はしていません。 そのうち何か良い題材を見つけて性能評価をしてみたいと思います。

2016年7月13日水曜日

代数的ブロック化多色順序付け法のオープンソース実装を公開しました!

前回の投稿から結構間が空いてしまいました。

今日は、Githubにて、新しいプロジェクトAlgebraicBlockMulticolorOrderingを公開しましたので、そのお知らせと解説記事です。

概要

今回公開した"AlgebraicBlockMulticolorOrdering"は、岩下・高橋・中島(2009)『代数ブロック化多色順序付け法による並列化ICCGソルバの性能評価』、情報処理学会研究報告、Vol.2009-HPC-121 No.11による「代数的ブロック化多色順序付け法」を実装したものになります。 詳しくは後述しますが簡単に言うと、この手法を使うと、ガウスザイデル法を並列化する時に元の格子の幾何形状を知らなくても収束性劣化を抑えることができるようになります。

下の図は、代数的ブロック化多色順序付け法を導入した時の効果を表したものです。 横軸が反復回数、縦軸がその時の残差(0だと完全に理論解と一致)を表しています。ですので、反復するほど残差が減少していきます。 ここで、緑色のSymGSというのは「対称ガウスザイデル法」というガウスザイデル法の1種であり、理想状態だと思ってください。 これを、通常の多色順序付け法(MulticolorOrdering)によって並列化した場合が茶色の線になります。ご覧のように同じ反復回数でも残差が大きく、収束性を劣化させている事がわかります。 一方で、ブロック化多色順序付け法(BlockMulticolorOrdering)を使った暗緑色の線では、収束性の劣化が抑えられている事が分かると思います。

このような効果をもたらすものが、代数的ブロック化多色順序付け法で、今回はこの手法を実装しオープンソースとして公開しました。 以下、この手法について「ガウスザイデル法をどのように並列化するか」という観点から、簡単ですが解説したいと思います。

ガウスザイデル法

ガウスザイデル法そのものについては、理学部・工学部であれば大学学部で習うと思いますが、ざっとおさらいしておきます。 「もちろん常識だよ!」という方は次の章まで飛ばしてもらって構いません。

ガウスザイデル法とは、線形方程式系の求解法の一種です。線形方程式系とは、連立一次方程式とも呼ばれている、中学校でも習う \[ \begin{align*} 3x + 2y & = 1 \\ 5x + 8y & = 9 \end{align*} \] みたいなものです。ここでは一般化して、 \[ \sum_j a_{ij} x_j = b_i \] と書くことにします。ここで、\(a_{ij}\)は↑の例だと3,2,5,8の係数のことで、\(x_j\)は\(x,y\)の未知変数、\(b_i\)は1,9の右辺定数に該当します。 \(\sum_j\)は\(j\)について総和するものです。

今、未知変数がn個あるとすると、方程式もn個あるはずな(というか無いと解けない)ので、iとjは0からn-1までの範囲の整数をとります。n=2の場合に書き下すと \[ \begin{align*} a_{00} x_0 + a_{01} x_1 & = b_0 \\ a_{10} x_0 + a_{11} x_1 & = b_1 \end{align*} \] となります。式が2つあり、それぞれ上から順に式0と式1と呼ぶことにします。

さて、この線形方程式系を解く方法は世の中にたくさんありますが、ここでは何らかの反復式を繰り返すことで解を収束させることを目的にした反復法を考えます。 そのためには、↑の\(a_{ij} x_j = b_i\)を\(x' = f(x)\)みたいな形にしなければなりません。

簡単にぱっと思いつくのは、式iを変数\(x_i\)の反復式に変形することです。つまり、 \[ a_{ii} x_i + \sum_{j \ne i} a_{ij} x_j = b_i \] と、\(x_i\)の部分だけ取り出してやって、あとは変形して \[ x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j \ne i} a_{ij} x_j \right) \] という式を出します。 これを反復式にするためには、右辺を既知の値とすれば良いので、k番目の反復時の\(x_i\)を\(x_{i, k}\)と表記すると \[ x_{i, k} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \ne i} a_{ij} x_{j, k-1} \right) \] という簡単な反復式ができます。これを繰り返せばなんとなく収束しそうですね。コードで書くと以下のようになります。

void Jacobi(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(n)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x_prev[j];
            }
        }
        x_next[i] = x_i / a[i][i];
    }
}

この方法はヤコビ法と呼ばれていて、最も簡単で、かつ各式を並列で計算できるという長所があります。 ただしこのヤコビ法には欠点があって、なにしろ収束がとても遅いのです。下手するとどこかで飽和して収束しないということにもなります。

これを改良するための方法が、ガウスザイデル法です。 方法は簡単です。今、i番目の変数\(x_i\)の計算時には、右辺は全てひとつ前の反復の値を使っています。 しかし、逐次実行していれば、i番目の式を計算している時には既に0からi-1の今の反復した後の値が分かっているはずですので、それを使うのです。つまり \[ x_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j < i} a_{ij} x_j - \sum_{j > i} a_{ij} x_j\right) \] とします。コードで書くと以下のようになります。

void GaussSeidel(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < i; j++)
        {
            x_i -= a[i][j] * x_next[j];
        }
        for(auto j = i + 1; j < n; j++)
        {
            x_i -= a[i][j] * x_prev[j];
        }
        x_next[i] = x_i / a[i][i];
    }
}

あるいは単純にxを1つにまとめて以下のように書くほうが簡潔で分かりやすいかもしれません。

void GaussSeidel(const double** a, double* x, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
}

この結果、下の図のように、収束性が格段に上がります。

ヤコビ法(Jacobi-青色)とガウスザイデル法(GaussSeidel-赤色)の比較

このように、ガウスザイデル法の方が収束性が良い一方、欠点もあります。 それは、先に述べた通り、この方法があくまで逐次実行することを前提としているという問題です。つまり、このままだと並列化できません

ということで、太古の昔は、並列性を取るならヤコビ法、収束性を取るならガウスザイデル法、という時期もあったようです。 しかしなんとかしてガウスザイデル法も高速化したい・・・そのための方法が、次に説明する多色順序付け法です。

多色順序付け法

多色順序付け法は、Multi-color OrderingやMCと呼ばれるものです。

多色順序付け法の説明に入る前に。まず、先の線形方程式系 \[ \sum_j a_{ij} x_j = b_i \] は、 \[A \boldsymbol{x} = \boldsymbol{b} \] と、係数行列\(A\)、未知変数ベクトル\(\boldsymbol{x}\)、右辺定数ベクトル\(\boldsymbol{b}\)で表せることは簡単に分かると思います。

ここで、線形方程式系が現実世界でどんな状況で出てくるか考えてみると、大体が\(A\)は疎行列になることがほとんどです。疎行列とは、\(A\)のうち殆どの要素が0であるような行列です。 どういうことかというと、例えば、熱伝導方程式(あるいはラプラス方程式) \[ \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} + \frac{\partial^2 \phi}{\partial z^2} = 0 \] を解きたい時、2階偏微分が \[\begin{align*} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i,j,k} & = \frac{\phi_{i-1, j, k} - 2 \phi_{i, j, k} + \phi_{i+1, j, k}}{{\Delta x}^2} \\ \left. \frac{\partial^2 \phi}{\partial y^2} \right|_{i,j,k} & = \frac{\phi_{i, j-1, k} - 2 \phi_{i, j, k} + \phi_{i, j+1, k}}{{\Delta y}^2} \\ \left. \frac{\partial^2 \phi}{\partial z^2} \right|_{i,j,k} & = \frac{\phi_{i, j, k-1} - 2 \phi_{i, j, k} + \phi_{i, j, k+1}}{{\Delta z}^2} \\ \end{align*}\] という形で中心差分で離散化できるので、3次元空間を正方形で区切った時に各頂点の値を求めるのには、上下左右奥手前の6点とだけつながっているモデルを作ることができます。 そのような時にそこから生成される線形方程式系は、1つの式が1つの頂点に対応しており、つまり行列で言えば、1行が1点に対応します。 とすると、1点はたかだか6点としかつながっていないことから、1行は6列だけ非ゼロで、あとは全てゼロとなります。

熱伝導方程式とはちょっと違いますが、平面で周囲8点(上下左右斜め)から影響を受けるという、もうちょっとまともに離散化したモデルだと行列は以下のような感じに見えます。

この行列はサイズが64x64の正方行列で、対角項が8(紫色)、非対角項には8点だけ-1が入っており(赤色)、他の要素は0(白色)になります。

さて、ということは、ガウスザイデル法の計算方法を再掲しますが

void GaussSeidel(double* x_next, const double** a, const double* x_prev, const double* b, const std::size_t n)
{
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < i; j++)
        {
            x_i -= a[i][j] * x_next[j];
        }
        for(auto j = i + 1; j < n; j++)
        {
            x_i -= a[i][j] * x_prev[j];
        }
        x_next[i] = x_i / a[i][i];
    }
}

で、ある行x[i]を計算するにあたっては、実はたかだか8個のx[j]からしか影響を受けません。 例えば、0点(行)目の計算では、1,8,9点(列)目からしか影響を受けないことが分かります。 また、2点(行)目は、1,3,9,10,11点(列)目からしか影響を受けません。 ということは、0行目と2行目はお互いに影響しておらず、実は並列で計算しても問題ないということが分かります。

このようにお互いに影響しない点同士を集めてグループにして、そのグループの中では並列で実行する方法を多色順序付け法と呼びます。 順序付け(ordering)と呼んでいるのは、計算順を上から順ではなくて、お互いに影響しない点のグループ順にする、という意味です。

多色(multicolor)の方は、実は歴史的経緯で、元々は、平面で上下左右の4点からしか影響を受けないような問題について考えた時に、 以下の図のように、黒色の点を計算する時には赤色からのみ影響を受け、赤色は黒色からのみ影響を受ける、と点を赤と黒の2グループに分けたことに起因します(この場合だけを赤黒法、Red-Black GaussSeidelとも呼ばれます)

この多色順序付けですが、単に計算順序を変えるだけだと、飛び飛びの行を計算することになり、キャッシュ効率等が悪くなってしまいます。 そこで、実際には、計算順序に沿って行列・ベクトルそのものを並び変えることが多いです。

下の図は、先ほどの行列を多色順序付けして並び替えたものです。

若干分かりづらいのですが、この行列は4色(グループ)に分けることができました。元々のサイズが64だったので、16行(点)ずつに分かれたことになります。 実際、先頭の色である0-15点(行)目は、同じ色の0-15点(列)目との係数はすべて0になっていることが分かります。次の16-31行目も同様です。 なので、この行列は、16行ずつなら並列化することが可能、ということになります。 16という絶対数だけ聞くと少なく聞こえるかもしれません。 しかし、問題サイズが大きくなっても、点と点の接続状況が変わらないかぎり色の数も同じなので、例えば32万点だった場合、8万点ずつ並列で計算できることになるので、計算時間がつらくなる問題サイズが大きいものでは十分に効果を発揮します。

多色順序付け法の欠点

ということで、ガウスザイデル法を多色順序付けを使うことで並列化することができました。めでたしめでたし。 と言いたいところなのですが、実は多色順序付けすると、元の逐次実行に比べて欠点があったりします。 それはベクトルデータのランダムアクセス化と、収束性の悪化です。

ベクトルへのランダムアクセス

元の行列と多色順序付けした後の行列の様子を並べて再掲します。

左:元の行列、右:多色順序付けして並び替えた行列

ここで、ガウスザイデル法は、基本的には疎行列・ベクトル積と同じく、係数行列\(a_{ij}\)に入出力ベクトル\(x_j\)をかける演算でした。 つまり、各行は、その行の非ゼロ(図でいうと色の付いている)部分に該当するベクトルxの値をアクセスします。

ということは、元の行列だと、例えば、0行目は、自分自身(0)と1,8,9列目のベクトルxの値が必要で、1行目は0,1,2,8,9,10列目が必要です。 という感じで見ていくと、元の行列のままだと、i行目の計算をする時には、基本的には中央のi-1,i,i+1列目と、そこから6つずつ左右に離れた3列ずつ、という3つの塊でアクセスしていることが分かります。 このように固まっている上、上から少しずつずれていっているので、隣の行を計算するときには既にキャッシュに乗っている可能性が高いです。 つまり、キャッシュには割と当たりやすい形をしています。

一方で、多色順序付けをして並び替えた方は、例えば0行目は0,16,32,48列目となり飛び飛びの場所のベクトルの値を必要とします。 1行目も同様に1,16,17,33,48,49と飛び飛びです。 これは当然といえば当然で、なぜなら多色順序付けは並列化のために格子の上で隣り合った(接続している)点を別のグループに分割するという手法だからです。 なので、今回の例のように4色に分かれていれば、基本的に必ず他の色のどこかにはアクセスする必要があるので、離れた4箇所にアクセスする必要が出てきます。 一般的に、n色に分けられた場合はベクトルはn箇所に分散します。 そのせいでキャッシュ効率が下がり、メモリアクセス速度がとても遅いなど、場合によっては並列化すると逆に遅くなる、ということも起きたりします。

というわけで、多色順序付けは並列化できる一方で、原理的にキャッシュ効率の悪化を避けられない、という問題があります。

収束性の悪化

キャッシュ効率が落ちるのは問題ですが、どれほど落ちるのかはやってみないと分からない上に、その対象とする問題や実効環境によっては別にあまり気にならなかったりします。

一方、こちらの収束性の悪化は結構悲惨です。 収束性の悪化も、最終的には、やはり元の問題(行列)の性質によるのですが、基本的にはどんな問題でも少しは悪化するはずです。

収束性が顕著に悪化する例として、ここでは対称ガウスザイデル法を例に挙げます。 通常のガウスザイデル法を見て気づいた方もいると思いますが、基本的には後に計算した値の方が前に計算した値より改善された(より解に近い)値を得ることができます。 ということは、0行目の結果とn-1行目の結果では、明らかに後者だけがとても解に近く、前者は離れたところにいるという非対称性を生むことになります。

非対称になると困る場合はいくつかあり、これを解説するとそれだけで記事1つできてしまうので詳しくは述べませんが。 ただ、多くの場合は、ガウスザイデル法をなんらかの他の反復法の前処理として使う場合に、前処理後に対称性が崩れてほしくないということに起因して対称性を残したいことが多いです (※たまに、対称ガウスザイデル法はガウスザイデル法より収束性が高いと勘違いしている人も多いですが、それは嘘です。 対称ガウスザイデル法は↑のように2回ループが回るので1回のガウスザイデル法より収束性が高いように見えがちですが、実は素直に2回ガウスザイデル法を回すほうが収束性が高いです)

このような非対称性を改善するものが以下の対称ガウスザイデル法です。

void SymmetricGaussSeidel(const double** a, double* x, const double* b, const std::size_t n)
{
    // forward sweep
    for(auto i = decltype(n)(0); i < n; i++)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
    // backward sweep
    for(auto i = static_cast<std::make_signed_t<decltype(n)>>(n)-1; i > 0; i--)
    {
        double x_i = b[i];
        for(auto j = decltype(i)(0); j < n; j++)
        {
            if(j != i)
            {
                x_i -= a[i][j] * x[j];
            }
        }
        x[i] = x_i / a[i][i];
    }
}

この対称ガウスザイデル法では、0行目からn-1行目まで順に計算した後、逆にn-1行目から0行目に向かって逆順にも計算する手法です。 このように逆回しも実行することで、どの値もちょうど同じぐらい改善されており、結果が対称になります。

と、対称ガウスザイデル法の話で少し脱線してしまいましたが、では実際に多色順序付けをするとどれほどガウスザイデル法の収束性が劣化するのかを示したのが下の図です。

まず、なにもしないガウスザイデル法(GaussSeidel、赤色)と多色順序付け(MulicoloOrdering-GS、黄緑色)はそれほど収束性に差がありません。 一方で、なにもしない対称ガウスザイデル法(SymGS、緑色)に対して多色順序付けした場合(MulicoloOrdering-SymGS、暗赤色)は大きく収束性が劣化して線が随分上に行ってしまっているのが分かります。 なお、ついでなのでガウスザイデル法を素直に2回まわした結果も乗せました(GaussSeidel x2、黄色)。 ご覧のとおり、SymGSよりこちらのほうが収束性が高いことが分かります。

今回の例ではこういう結果になりましたが、先述の通り実際には、多色順序付けすると普通のガウスザイデル法だけでも収束性がもっと劣化する場合もあります。 逆に、対称ガウスザイデル法でもあまり劣化しないこともあります。 が、いずれにせよ多色順序付けすると多かれ少なかれ収束性が劣化するという性質があることは間違いありません。

それは、多色順序付けで収束性が劣化する原因が、主にはincompatible nodeの増加によるものだからです。 incompatible nodeとは、「周囲の点/行がひとつも改善されていない状態から始める点/行」のことです。 詳細な説明等は中島先生の資料等に譲りますが、直感的には、なるべく周囲に改善された(解に近い)値があったほうが自分自身の解はより良いものになるので、その数が少ないと収束性が劣化する、ということです。

なぜincompatible nodeが増えるのかというのは、これも多色順序付けが何をしているのか考えれば当然で、同じグループの中の行は互いに独立して計算ができるようになっているので、 最初に計算する色(グループ)に属する行は全部周囲の点はまだ一度も計算されていないことになり、incompatible nodeに該当してしまうからです。 対して、何も並び替えず逐次で実行する場合は、基本的には最初の1行/点だけで、後続の処理は基本的には全部、周囲のどこかは既に改善された点になっているはずです。 つまり、n次元正方行列がc色に多色順序付けされたとすると、incompatible nodeの数は必ず1点からn/c点に増加します。

ということで、多色順序付けをすると、原理的に収束性を劣化させることにつながる(かもしれない)のです。 せっかく多色順序付けをして並列化できるようにしたのに、収束までに必要な反復回数が増えてしまってはあまり効果がありません。

ブロック化多色順序付け法

さて、通常の多色順序付け法ではこのように原理的に避けようがない問題が2点ありました。さてどうしましょうか。これを避ける方法がブロック化です。

ブロック化の説明のために、先に紹介した2色で塗り分けられる4点ステンシルの以下の例を挙げて説明します。

この場合、黒→赤の順で計算が実行されるとすると、この状態だと黒色が全てincompatible nodeになり合計で8個あります。 また、黒色が必要とする周囲の赤色の点は全て遠いところにあるのでキャッシュ効率も下がっているという先ほどの欠点があります。

しかし、これを、2x2の点でブロック化して見てください。つまり行番号のままで言うと

  • 左下のブロック#0:0,8,10,2
  • 右下のブロック#1:1,9,11,3
  • 左上のブロック#2:4,12,14,16
  • 右上のブロック#3:5,13,15,7

の4つのブロックに分けます(※ブロック≠色のグループです、気をつけてください)。 そうすると、ブロック単位でみると、左下のブロック#0は右上のブロック#3とは繋がっていないので独立して計算できることが分かります。また、右下のブロック#1と左上のブロック#2は独立しています。

と、いうことは、ブロックに分けた状態でも独立しているブロック同士は並列で計算できます(※ブロックの中の点同士は接続されているので、逐次実行する必要があります)。 もうちょっと言うと、ブロックに分けた後に、ブロック間の接続状況で多色順序付けすれば、ブロック同士は並列計算できるようになります。これが、ブロック化多色順序付けです。

ブロック化の効果

ブロック化すると、先に上げた多色順序付けの問題が緩和されます。

先の例だと、まずincompatible nodeはブロック化することによって8から4に半減します。 なぜなら、ブロックの中は逐次実行されるので、基本的にはincompatible nodeは「最初に計算される色に属するブロックの中の最初の点」だけになるからです。 例えばどのブロックも左下から計算するのであれば、incompatible nodeは、以下の4点だけになります。

  • ブロック#0の点0
  • ブロック#1の点1
  • ブロック#2の点4
  • ブロック#3の点5

実際に収束性が改善されている様子は、先に通常の多色順序付け法の収束性を論じていた時と同じ対象にブロック化を適用して実験してみた以下の結果を見ればよく分かると思います。

ブロック化していない普通の多色順序付け(MulticolorOrdering-SymGS、暗赤色)に比べて、ブロック化多色順序付け(BlockMulticolorOrdering-SymGS、水色・暗緑色)は、 理想的な逐次板対称ガウスザイデル法(SymGS、明緑色)にぐっと近くなっている様子がわかります。 また、2x2の4点ブロック(水色)にするより、ブロックのサイズ(1ブロックに含まれる点の数)を4x4に増やした方(暗緑色)が、より効果的であることが分かります。

ブロック化すると行列はどうなっているかというとこんな感じ↓になっています。

この例は2x2(つまり4点)でブロック化したもので、ブロックごと4色に分割されています。 この行列の様子を見ると分かる通り、上から順に実行したとして、周囲にまったく改善された点が含まれていない点(incompatible node)に該当するのは、0,4,8,12行目だけです。 なにもしない多色順序付けだと最初の色(0-31行目)が全てincompatible nodeで32点あるのに対して、4点まで減らすことに成功しています。 そしてまた、ブロックの中は接続されている上、他のブロックの点を必要としないor必要とする色の数が少ない点があるので、キャッシュ効率も上げられていることが分かります。

一般的に、n次元正方行列をブロックサイズbにブロック化してc色に分けられた場合、incompatible nodeはn/(bc)にできます。なので、可能な限りブロックサイズは大きくしたほうが得策ですね!

ブロック化の欠点

と、ブロックサイズは大きくしたほうが良いと言いましたが、半分嘘です。

勘の良い方は気づいていると思うのですが、ブロック化したことによって、並列実行できる数(並列度)が減っています。 今回の例だと、ブロック化しなければ1色の中の全ての行、つまり16行を並列実行でき、前回の記事に書いた通り、一般的に一度に並列実行できる数はn/cになります。 一方、ブロック化するとブロック間でしか並列実行できないので、4ブロックずつです。一般的に並列実行できる数はn/(bc)です。 そう、incompatible nodeの数と並列実行できる数は基本的には一致します。 極端な話を言えば、n次元正方行列をブロックサイズnにするとincompatible nodeは1にできますが、全く並列化出来ません。ブロックが正方行列全体になってただのガウスザイデル法になってますので。

ただ一方で、そもそも(効果的に)並列実行できる数は計算するハードウェアによって上限が決まっています。普通のCPUだと多くても16か32ですし、GPUでもせいぜい1kか多くて10k個オーダーです。 なので、並列計算が必要になるような問題サイズが十分に大きい場合は、ブロック化しても速度はあまり低下しないと考えられます。 理想的には、n次元正方行列がc色に分けられるのであれば、最適な並列数がpのハードウェア上で計算する時はn/(bc)=pつまり、ブロックサイズbはn/(cp)にすると良いです。

とは言うものの、ブロックサイズを大きくし過ぎると、隣の並列コアと処理する場所が遠くなるので、必ずしもブロックサイズをこの理想数にしたほうが良いというわけでもないようです。 最終的にはキャッシュの大きさ等と相談して決める事になると思いますが、キャッシュ効率はともかく、収束性はどれだけブロックを大きくしても逐次板に近づくだけなので一定以上は上がりません。 経験的には、ブロックサイズは大きくて1行あたりの非ゼロ要素数の3割ぐらいで十分なようです。

代数的ブロック化多色順序付け法

さて、ブロック化できましためでたしめでたし・・・と行きたいところなのですが、 ↑の方法は、格子の状態を見て説明をした通り、基本的は行列を生成する元の問題の形状に依存してブロックに分類します。 そのため、幾何的ブロック化多色順序付け(Geometic Block Multicolor Ordering, GBMC)と呼ばれます。

でも、元の問題なんてわからないことも多いですし、元の問題によってブロック化の方法を変えるのって面倒ですよね。 そこで、行列の接続状況からだけでブロック化する方法があります。 それが代数的ブロック化多色順序付け(Algebraic Block Multicolor Ordering, ABMC)です。 今回公開したものは、このABMCを実装したものになります。

ABMCの方法については、先述の岩下ら(2009)の論文に詳しく書かれています。 英語でよければ、同じ著者による、IEEE Trans.の論文スライド資料もあります(発表年順からすると、IEEE Trans.の方が元のようです)。 詳しくはこれらの論文等を読んでもらうとして、おおざっぱな方法は以下のようになっています。

  1. まだどのブロックにも入れていない行(点)iを探す
  2. 点iを、新しいブロックbに入れる。
  3. 点iに接続している列(点)のうち、まだどのブロックにも入れていない列(点)jを探す
  4. 点jをブロックbに入れる
  5. ブロックbが一杯になったら最初から繰り返す
  6. まだ空きがあるなら、今度は点jに接続している点kを探す
  7. 点kを新しいブロックbに入れる
  8. ブロックbが一杯になったら最初から繰り返す
  9. 以下、今度は点kに接続している列(点)を探してまたbに入れて・・・を繰り返す
  10. 全点がブロックに入ったら、今度はブロック間の接続状況を確認して、普通の多色順序付けをすれば、完成

今回公開したコードでは、Abmc.hppの321行目からの関数がこの部分に該当します。 文章で説明されるよりコードを見たほうが分かりやすい方はぜひそちらをご覧ください。

まとめ

ということで、ガウスザイデル法の復習から、代数的ブロック化多色順序付け(ABMC)法までを簡単ですが解説しました。

この法に限らず、論文に手法自体は載っていても、実際の実装は公開されていることはほぼ稀で、仕方ないので自分で実装する時に苦労することがよくあります(特に数値計算分野)。 今回このコードを公開することで、後に続く方々の参考になれば、とても嬉しく思います。

なお、今回のコードは、これ自体は特に高速化等をあまり考えて作ってありません。 ひとつは先述のように「他の方の参考に」という意図を含めて分かりやすさを重視しているのと、またそもそも行列の前処理部分は実際の計算部分より小さく無視して良いことが多いから、という理由からです。 ただ、とても挑戦しがいのある問題だとは思いますので、またそのうち機会があればやってみた結果をご報告させていただくことができるかもしれません(アイディアはあるので)。 もし、並列化実装してみたとか、あるいはそれ以外でもバグ等見つけられた方は、遠慮無くissueやpull requestを投げてください。 お待ちしています。

あと最後に一応宣伝です。 今回のような「アルゴリズムそのものを見直す」という部分も、フィックスターズの得意とする「高速化」に含まれており、アルゴリズムからハードウェア特化したものまで、総合的なな解決方法を提供しています。 もし今現在何か「これを速くしてみたいんだけど」とお困りのプロジェクトがありましたら、ぜひお問い合わせください。 また、そういった仕事に興味がある方は、ぜひ一緒に働いてみませんか?

2016年3月13日日曜日

C++のムーブと完全転送を知る

社内勉強会、今期(と言ってももうあと1ヶ月もないですが)は、数理最適化勉強会と、Effective Modern C++輪読会をしています。この記事は、後者のEffective Modern C++輪読会で、『Effective Modern C++』5章の一部を輪読した時の資料を流用したものです。

C++11と言えば、昔のC++(03)から色々あって多くの機能が追加されとても便利になったバージョンです。さらに、C++14は11では間に合わなかった・忘れていた色々な便利なものを補填したもので、Effective Modern C++輪読会は、『Effective Modern C++』を教科書にしながらこのC++11/14について学ぶ会になっています。

時は既に2016年、gccもclangもMSVC++も概ねC++14が使えるようになっており、もはやC++14が使えないコンパイラにはC++コンパイラを名乗る資格はないと言っても過言ではない時代です(※個人的な意見であり弊社エンジニアの総意ではありません)。C++14が現代機械だとすると、C++11は戦後直後ぐらいの道具、C++03あたりは江戸時代ぐらいの遺物という感覚でしょうか(※あくまで個人的な感覚です)。

そんな状況で、Effective Modern C++輪読会と言うとちょっと今更感はあるのですが、こういうものはなんとなくの理解で進んでしまうことが多いので、この勉強会では学ぶことはけっこうありました。特に、今回のstd::movestd::forwardは、C++11で導入された機能の中で便利だと聞くけどいまいちよく分からないものNo1(※個人比)です。既に世の中にはもっと優れた記事がたくさんあるのは確かですが、色々な観点からの説明があるに越したことはないと思っているので、せっかくなので記事にして残しておきたいと思います。

一方で、C++といえばマサカリを投げ合う文化で有名ですが(?)、この記事も、大きな間違いはないとは思ってますが、マサカリを投げる隙がまったくない自信はありません・・・(先述の通り難しいところだと思いますし)。もし何か誤りなど見つけられた方は、記事にコメントいただくか、Twitter:@LWisteriaまで直接ご連絡ください。

用語の復習

ムーブと転送の説明に入る前に、用語の軽い復習をしておきます。

右辺値と左辺値

式には右辺値(rvalue)と左辺値(lvalue)があります。厳密な定義は複雑で難しいですが、概ね「アドレスを取得できる式=左辺値、できない=右辺値」と覚えておけば良いです。

int z = f(x, g(y));

ここで、

  • x, y, z:左辺値
  • g(y), f(x, g(y)):右辺値

です。

あるいは、これも厳密ではないですが、左辺値は名前のついた変数、右辺値は一時オブジェクトのように覚えてもいいです。

左辺値・右辺値それぞれの参照は左辺値参照(lvalue reference)と右辺値参照(rvalue reference)で、*通常は*、前者がT&で後者がT&&と書かれます。

int& y = x; // 左辺値参照
int&& y = g(x); // 右辺値参照

ここであえて*通常は*と言ったのは、T&&で宣言しても左辺値になることがあるからです。

void f(int&& v)
{
    // vはT&&(右辺値参照)で宣言しているのに、アドレス取得できるので、実はvは左辺値です!
    // ∵すべての関数の引数は左辺値だから
    std::cout << &v << std::endl;
}

このようになる詳しい仕組みは後述しますが、この先、説明中にT&&が出てきても右辺値とは限らないことがあることは念頭においておいてください(逆に、T&&が毎回出てきた時に、毎回、それが右辺値か左辺値かを考えていると、理解が進むかもしれません)。

ムーブとコピー

ムーブ(move)は、コピー(copy)の対になる概念で、コピーが複製なのに対して、ムーブは置換(複製&破棄)を意味します。 RPGツクールで物体の移動をやりたい時に「移動」というコマンドがなくて「複製してから消せばいい」と書いてあったあれがムーブです。

このムーブを実現するコンストラクタおよび代入演算子をムーブコンストラクタムーブ代入演算子と呼びます。

{
    const int y = x;
    const int z = x;
    int ret = 0;
    ret = y+z;
}

↑の時、代入演算=のうち、y=xの=はコピー代入である必要があります。一方、z=x, ret=0, ret=y+zの=は、(コピー代入でもいいですが)xy+zの結果を捨ててもいい(後で使わない)のでムーブ代入にしてもよいはずです。

このように、同じ「代入」という演算は、コピーの意味とムーブの意味があります。「ムーブの意味にすること」をムーブセマンティクス(で解釈する|をとる)、と言います(※実はあんまり厳密な説明はできてないので、yohさんの『本当は怖くないムーブセマンティクス』(C++アドベントカレンダー2012 15日目)あたりを参考にしてください)。

(完全)転送

転送(forward)とは、受け取った引数をそのまま別の関数に受け渡すことです。

void f(int x)
{
    g(x);
}

そして完全転送(perferct forward)とは、仮引数の値だけでなくその他の情報(型)も含めて完全に実引数に渡すことです。↑の例だと、f(z);と呼び出した時に、zの情報をそのままg()に伝え、g(z);と呼び出したことと等価になるような転送の場合を、完全転送と呼びます。

ですので、この例では完全転送を実現できていません。もしint& z;と宣言されていたら、f(z)g((int)z)と勝手に参照&を外されていることになるからです。

ムーブ

std::moveとは

C++でムーブを実現するための機能がstd::moveです。しかし、その名前に反してstd::moveはムーブしませんstd::move(x)はxをムーブできるようにするだけです。

この「ムーブできる」とはどういう状態かというと、先に上げた例だと

{
    const int y = x; // ムーブできない
    const int z = x; // ムーブできる
    int ret = 0;     // ムーブできる
    ret = y+z;       // ムーブできる
}

となっており、ここで右辺に着目すると、左辺値xはムーブできない場合もあるが、右辺値0y+zは必ずムーブできています。 それもそうで、ムーブできるのは、右辺を廃棄して良い時なので、右辺値(≒一時オブジェクト)なら必ずムーブ可能ということです。

つまり、std::move右辺値にキャストすれば実現できます。実際に使ってみるとこんな感じです。

    f(x+1);          // x+1は右辺値
    f(x);            // xは左辺値
    f(std::move(x)); // moveするとxが左辺値になる

ここで、std::moveは実際にはムーブせずにムーブできるようにする(キャストする)だけと言っているのは、つまりstd::moveしてもムーブにならないことがあるからです。

どういう時かと言うと、const値をstd::moveした時です。

static void Foo(const Pointer p)
{
    Pointer q(std::move(p)); // moveしたのに、コピーコンストラクタが呼ばれる!
}

なぜこんなことが起きるかと言うと、ほぼ全てのムーブコンストラクタ・代入演算子は、非const参照を引数とするからです。 もう少し詳しく言うとムーブとは「代入元を破棄する」処理のことなので、引数(代入元)を操作(破棄)するため、非const参照を必要とします。

    // コピーコンストラクタ
    Pointer(const Pointer& f)
    {
        this->ptr = f.ptr;
    }
    
    // ムーブコンストラクタ
    Pointer(Pointer&& f)
    {
        this->ptr = f.ptr;
        
        // ムーブしたあとは廃棄する
        f.ptr = nullptr;
    }

なので、私のような手が勝手にconstをつけるような人は特にconst値はムーブできないよと手に教えこむ必要がありますので、注意しましょう。

なお、「ほぼ」と書いたのは、規格上は別にconst右辺値参照を引数にすることもできるからですが、それはもはやただのコピーなので意味がないです。 (※復習:左辺値参照をとる関数に、右辺値を渡しても問題ない。逆に、右辺値参照をとる関数に、左辺値は渡せない。なので、大体の場合、const右辺値参照を取る関数は不便なだけで、const左辺値参照を取る関数があれば済む)

std::moveの正体

先述の通り、std::moveの正体はただ右辺値にキャストするだけの関数です。 ではどうやって右辺値にキャストしているかというと、こんな感じになります。

template<typename T>
decltype(auto) move(T&& param)
{
    return static_cast<std::remove_reference_t<T>&&>(param);
}

やっていることは

  1. 関数の戻り値は右辺値になるので、同じ値を返すだけの関数にしたい
  2. 戻り値を参照にしないとコピー(あるいはムーブ)が走るので、autoではなくdecltype(auto)
  3. 同様に、引数も参照にしたいが、T&だと左辺値しか受け取れないので、T&&にする(※T&&は後述の通り右辺値参照とは限らない)。
  4. 最終的に、右辺値参照を返したいが、もしTが左辺値参照だとdecltype(param)(つまりT&&)も左辺値参照になってしまうので、確実に右辺値参照(&&)にキャスト

という流れです。

と、こんな感じで「ムーブできるようにする」意味からstd::moveという名前になりましたが、やってることは「右辺値にキャストする」なのでstd::rvalue_castみたいな名前も提案されていたようです。 個人的にはstd::movableとかが良かったんじゃないかなぁと思ってます。

ムーブできないこともある

このように、C++11以降はムーブが使えるようになったので、ムーブを使えば速くなる!と言われますが、そうじゃない場合もあるので過度に期待し過ぎないほうがよいと言われています。そうじゃない場合というのは、ムーブできない場合もあるからです。

まずそもそも、ムーブ演算(コンストラクタ、代入)は、明示的に実装してあるか、そうでないなら

  • コピー演算が明示的に実装されていない
  • デストラクタが明示的に実装されていない
  • 全てのフィールドはムーブ演算可能

の全ての条件が満たされた時にしか暗黙的に生成されません。 つまり、ムーブしたくなるような複雑で大きいものでは、だいたいどれかの条件にひっかかるので、明示的に実装していないとムーブ演算なんて存在しないことがほとんどです。

また、もしムーブ演算があっても、ムーブがコピーより軽いかどうかは実装依存です。例えば、std::vectorは確かにムーブはポインタの付け替えで済むので高速ですが、std::arrayは中身を複製しないといけないのでムーブもコピーも同じぐらいです。

遅くなることはまぁまずないですが、実装依存なのでそういうムーブ演算を存在させることは可能ですね。 なので、ムーブを使えばコピーより必ず速くなるとは思わない方が良いです。

さらに、ムーブ演算が存在し、確実にムーブがコピーより軽い実装でも、例外安全性のstrong保証(途中で例外が発生しても必ず操作前に巻き戻って返せる性質)を必要とするアルゴリズムでは、例外を投げない(noexceptがついた)ムーブ演算でなければ、ムーブ演算を使えません。

というように、ムーブ演算が効果を発揮しない場合を列挙するとキリがありませんので、あんまり期待しても仕方ないよね、という話でした。

とは言うものの、とあるコンパイラではC++11にしただけで1.x倍になったという事例もあるらしいので、やっぱりムーブのおかげで速くなることも結構あります。 どっちやねん!って感じになりますが、「期末試験よくできたかもと思ったけど点数低かったら傷つくから、点数は低いものだと思っておこう」ぐらいの感じでいればよいのかなと思います。

完全転送

std::forwardとは

std::forwardは「完全転送した呼び出しにする」関数です。 完全転送とは、最初に述べた通り、仮引数の元の実引数の完全な型情報も他の関数の実引数に渡すことです。

ここで、fが多重定義されていて、そのfを呼び、かつそれ以外の共通の処理をまとめたテンプレート関数Fを考えましょう。

C++03以前の場合、まず、

template<typename T>
static void F(T v)
{
    // ここには何か共通の処理をして
    f(v); // オーバーロードされた関数も呼ぶ
}

と書いてしまうと思いますが、これは完全転送できていません

なぜなら、F(x)と呼んだ時にxvにコピーされるため、fに渡る引数は別の実体を指します(中身はコピーされているので同じに見えますが)。 それにそもそも、コピーコストが無視できないオブジェクトだと、コピーが走るのは避けたいです。

では、コピーされるのが問題だから参照にすればいいでしょうか?

template<typename T>
static void F(T& v)
{
    f(v);
}

いいえ、これはそもそも失敗です。なぜなら、右辺値を渡せないからです。

さて、そこで、C++11からは右辺値も参照が使えます!

template<typename T>
static void F(T&& v)
{
    f(v);
}

これで右辺値が渡せます。ただし、まだ不完全です。

なぜなら、復習のところでも書きましたが、仮引数は必ず左辺値なので、Fを経由すると右辺値参照を取るfを呼べないのです。よし分かったじゃあstd::moveすればいいだろう、と思いつくかもしれませんが、それはそれで今度は左辺値参照を取るfを呼べません。

つまり、完全転送を実現するには、Fの実引数が右辺値ならfの右辺値も右辺値に、左辺値なら左辺値にと適切にキャストする必要があります。 それを実現するのがstd::forwardです。

template<typename T>
static void F(T&& v)
{
    f(std::forward<T>(v));
}

これで完全転送できました。めでたしめでたし。

つまるところ、std::forward<T>(v)とは、vの

  • 参照元が右辺値なら、引数を右辺値に
  • 参照元が左辺値なら、引数を左辺値に

キャストする関数ということです。std::moveと同じく、実はstd::forwardもただのキャストであり、std::foward自体が転送するわけではないです。

std::forwardの仕組みと参照の圧縮

std::forwardを使った完全転送

template<typename T>
static void F(T&& v)
{
    f(std::forward<T>(v));
}

で、std::forwardTを受け取って、参照元が右辺値か左辺値かを判断していると書きました。つまり、「Tに参照元が右辺値か左辺値かの情報が入っている」ということです。

ただ、別に新しい型みたいなのが追加されるわけではないです。実は、T&&の型推論をする時には、左辺値を代入されたらTは参照、右辺値を代入された非参照にするような仕様になっています。 例挙げると、

template<typename T>
static void F(T&& v)
{
    f(std::forward<T>(v));
}

int x;
F(x); // xは左辺値なので、Tはint&
F(1); // 1は右辺値なので、Tはint

といった感じです。

ただ、実際にTがどっちなのかを直接観測するのは難しいのですが、TypeDisplayerのテクニックを使うと一応知ることができます。

なので、std::forward<T>

  • Tが参照の時、vの参照元は左辺値だったと判断して、左辺値参照を返す
  • Tが非参照の時、vの参照元は右辺値だったと判断して、右辺値参照を返す

という挙動をすれば良いわけです。

ではそんな挙動をどうやって実現しているのでしょうか。そこで参照の圧縮(reference collapsing)の出番です。

参照の圧縮ってなんのこっちゃと思うかもしれませんが、簡単に言うと、「実はコンパイル時には参照の参照があって、参照の参照を普通の参照に変換しているんだ!(なんだってー(AA略」みたいな話です。 参照の参照というのは、普通はコンパイルできない

int& && y = x;

みたいなやつのことです。こんなの出てこないだろうと思うかもしれませんが、例えばテンプレートとtypedef/エイリアスを使って

template<typename T>
struct L
{
    using Type = T&;
};

みたいなののTに参照を入れると、コンパイル時には登場します。例えば、Tに右辺値参照(int&&)を入れると、Typeは右辺値参照の左辺値参照(int&& &)になりますよね。

このようなコンパイル時に発生する二重の参照を単独の参照に変換することが参照の圧縮です。 参照には右辺値と左辺値があるので、パターンは4通りありますが、以下のような規則で変換されます。

便宜的に、二重の参照をT<1> <2>と書きます(例えばint& &&は<1>=&, <2>=&&)。この時の規則は以下の通り。
<1>/<2> 左辺値参照& 右辺値参照&&
左辺値参照& 左辺値参照 左辺値参照
右辺値参照&& 左辺値参照 右辺値参照

つまり例で言うと、int&& &&=int&&であり、その他は全てint&になるということです。 実際に確かめてみると、確かにそうなってるのがわかります。

こうしておくと、std::forward

  • T&&Tに非参照(int)や右辺値参照(int&&)を入れれば、右辺値参照が作れる
  • 左辺値参照(int&)を入れれば、左辺値参照を作れる

ようになります。

ということで、std::forwatd

template<typename T>
decltype(auto) forward(std::remove_reference_t<T>& param)
{
    return static_cast<T&&>(param);
}

と実装できます。

どういう動きをするかというと、int x;の場合、

forward(x)と左辺値を渡した場合
  1. Tは参照(int&)になる
  2. remove_referenceにより参照が外されて、非参照(int)になって、paramの型は参照(int&)に戻る。
  3. T&&int& &&と左辺値参照の右辺値参照になり、圧縮されて左辺値参照int&になる
  4. なので、全体としては、左辺値を受け取った時はその左辺値参照を返す
forward(1)あるいはstd::forward(std::move(x))と右辺値を渡した場合
  1. Tは非参照(int)になる
  2. 既に非参照なのでremove_referenceは何もせず非参照(int)のままで、paramの型は参照(int&)になる。
  3. T&&int&&と評価され右辺値参照になる
  4. なので、全体としては、右辺値を受け取った時はその右辺値参照を返す

という挙動を実現しています。

なお、参照の圧縮(つまり参照の参照)が起きる場合は実は4パターンしかありません。そのうち

  • typdefやエイリアス
  • テンプレートでのT&&での型推論

は既に説明しました。

3番目は、実はテンプレートではなくauto&&の型推論でも発生するということです。

int x;
auto&& y = x; // これは左辺値を代入するので、autoがint&になって、yはint& &&が圧縮されてint&になる
auto&& z = 1; // これは右辺値を代入するので、autoがintになって、yはint&&のまま

これも直接観測は難しいですが、TypeDisplayerで確かめると確かにそうなっています。

このように、テンプレートやauto型推論する場合には、参照の圧縮によって、&&と書いても右辺値ではなく左辺値参照になることがあります。 最初に「&&は右辺値参照とは限らない」と書いたとおりです。このような型推論を伴う&&で、実際には左辺値参照になるかもしれない参照は、ユニバーサル参照などと呼ばれています。

あと最後の4番目は、decltypeでも起きます。 decltypeは引数に変数名を与えると変数の型になり、式を与えるとその左辺値参照になります。なので、int x;の時

  • decltype(x)は、int
  • decltype((x))は、int&

になります。ということは、

template<typename T>
void f(T v)
{
    decltype((v)) a = v;
}

f<int&&>(1); // これは、Tがint&&なので、decltype((v))は、int&& &で、参照の圧縮がなされてaの型はint&となる

ということが起きます

完全転送できない場合とその対策

さて、先ほど以下の様な完全転送を実装しました。

template<typename T>
static void F(T&& v)
{
    f(std::forward<T>(v));
}

でもこれは実はちょっと不完全です。なぜなら、引数の数が違う多重定義に対応していないからです。

extern void f(int x);        //こっちはF経由で呼べるが
extern void f(int x, int y); //こっちはF経由では呼べない

template<typename T>
static void F(T&& v)
{
    f(std::forward<T>(v));
}

引数の数も違う場合に対応するために、可変長にしましょう

template<typename... Ts>
static void fwd(Ts&&... params)
{
    f(std::forward<Ts>(params)...);
}

これでどんな場合でも完全転送できるようになりましたね!!めでたしめでたし!

と言いたいところですが、実はこれでも、f(???)fwd(???)と???に同じものを書いた時に挙動が異なってしまうような完全転送できない場合が5パターン知られています

  • 参照がとれない場合
    • ビットフィールド
    • 定義のない静的(コンパイル時)定数フィールド
  • 型推論に失敗する・型推論結果が期待と異なる場合
    • NULL
    • 多重定義された関数・テンプレート関数
    • 波括弧の初期化子

でも大丈夫です。一応ちゃんと対策もあります。1つずつ見てみましょう。

まず最初。ここまで見てきたように、完全転送には元の引数に対する参照をとる必要があります(参照がないと、コピーしなきゃいけなくなるので)。 右辺値参照によって、ほぼすべての参照がとれるようになりましたが、実は参照がとれないパターンが2つあります。

参照が取れないパターン2つのうち、簡単な方は、ビットフィールドです。

struct Foo
{
    std::uint32_t a: 28,
                  b: 4;
};

これは、多くの場合、参照が取れません。

Foo f;
std::uint32_t& a  = f.a; // 左辺値参照もとれない
std::uint32_t&& b = f.b; // 右辺値参照もとれない

なぜか?参照はハード的にはポインタと同じなのでアドレスが必要ですが、例えば「4bit後ろのアドレス」って取れないからです。とれるのは、アドレスがとれるのはビットフィールドを汎整数型に変換した先だけです。

なので、逆に言えば、汎整数型に変換した先のアドレスで事足りるような場合、つまり

  • 非const参照
  • 元の構造体が右辺値

である時は値を書き換えないので、参照をとっても問題ないとみなされます。

    // const参照はOK
    const std::uint32_t& a = f.a; 
    const std::uint32_t& b = f.b;

    // 元の構造体が右辺値なら、非const参照も可
    std::uint32_t&& c = std::move(f).a;

どちらも、変換された汎整数型の実体(アドレス)を指す参照になります。

ということで、もう少し厳密に言うと、元の構造体が左辺値の場合にビットフィールドに対する非const参照の場合に参照が取れないので、この場合は完全転送もできないです。

    Foo foo;
    f  (foo.a); // OK
    fwd(foo.a); // NG

これを回避するには、明示的にコピーを作ればよいだけです

    Foo foo;
    const uint32_t a = foo.a;
    f  (a);
    fwd(a);
    foo.a = a;

参照が取れないパターンのもう片方は、定義のない静的(コンパイル時)定数フィールドです。

struct Foo
{
    static const std::size_t N = 10;
    static const std::string S;
};
// const std::size_t Foo::N = 10; // 定義がなくてもOK
const std::string Foo::S = std::string("foo"); // こっちは非汎整数型なので、定義も必要

のように、汎整数型の場合、定義がなくてもFoo::Nを使うことができます。なぜなら、定数伝搬(const propagation)するので、定義(実体)がなくても、コンパイラがFoo::Nを10に置き換えなければならないからです。 なのでこの状態で参照を取ろうとすると、実体がない=メモリが割り当てられていないのでアドレスがとれず、つまり参照がとれない(コンパイルはできるけどリンクエラー)になります。

これ割とよくあって、私がよくやるのが、コンパイル時定数(constexpr)を完全転送しようとしてリンクエラーです。「コンパイル時定数なのに定義が要るとかなんでやねん!」ってなるので気をつけましょう(※一応、言われるがまま定義をすれば動きます)。

ただし、コンパイラによっては定義がなくても動いてくれるらしくて、それも規格上は許容されているらしいです(でも少なくともgcc/clangでは動かないです)。

以上のような参照が取れない場合以外、つまり参照がとれたとしても完全転送できない場合があります。 今までずーっと見てきたように完全転送はうまく型推論を駆使することで実現していますが、逆に言うと、型推論ができなかったり、推論結果が違うと、完全転送できません。そんなパターンが3つ知られています。

一番簡単なのはNULLを使った場合です。 NULLは実体は(だいたいは)0なので、可能な限りポインタではなく汎整数型(通常はint)に推論されます。

なので、ポインタの代わりにNULLを完全転送しようとすると、コンパイルエラーになります

static void f(int* adr);
fwd(NULL); // NG
f  (NULL); // NG

解決策はnullptr使えです。

fwd(nullptr);

「今時まだNULLを使ってるなんておっさんみたい!」と現代っ子に言われないよう、NULLはもうやめましょう。

さて、次のパターンは、多重定義された関数・テンプレート関数です。 関数ポインタ型変数に関数を代入すると、関数への関数ポインタが取れるのは

void (*p)(const int) = f;

f(1);
p(1);

と、太古の昔から変わらない現象ですが、多重定義された関数を関数ポインタとして取る場合、受け取る側の型が合う方を代入してくれるのもご存知ですか?

static void f(const int val);
static void f(const char str[]);

void (*p)(const int val)    = f; // これは上のf
void (*q)(const char str[]) = f; // これは下のf

なんかいつもの型推論みたいに、右辺から左辺を推論するのではなく、左辺から右辺を推論しているので若干気持ち悪いんですが、まぁこれのおかげで多重定義されていても、型さえしっかりしてあれば関数ポインタを区別して取れます。

これは関数の引数にした時も同じです

static void g(const int val);
static void g(const char str[]);

static void f(void (*p)(const int val));

f(g); // fが受け取るのはintを引数に取る関数しかないので、このgは上側のgになる

が、これを完全転送しようとする失敗します

f(g); // OK
fwd(g); // NG

なぜなら、fwdはどんな型でも受け取る関数なので、いつものように代入先からgを推論することはできないからです。

回避策は、キャストするなどして代入元で型情報をつけてやればOKです。

fwd(g); // NG
fwd(static_cast<void(*)(const int)>(g)); // OK

なお、この問題はテンプレート関数でも同じです。

    template<typename T>
    static void g(const T val);

    f(g);
    fwd(g); // NG
    fwd(static_cast<void(*)(const int)>(g)); // OK
    fwd(g<int>); // OK

最後は波括弧初期化子です。これも、先の関数ポインタのように、代入先によって型を決定します

    std::vector<int> v = {0, 1, 2};
    auto&&           i = {4, 5, 6}; // initializer-list

ということで、完全転送しようとすると、関数ポインタと同じように、どの型になるべきか推論できないという問題が発生します。

    f({0, 1, 2}); // initializer-list
    fwd({4, 5, 6}); // NG

そして、この回避策も関数ポインタと同じで、ちゃんと明示的に型を指定すれば良いだけです。

    auto i = {4, 5, 6}; // こうすればinitializer-list
    fwd(i);
    fwd(std::vector<int>({4, 5, 6}));

まとめ

ということで、C++のムーブおよび転送についてざっと解説してみました。

実際使ってみないとややこしいだけに思えますが、特に高速かつ便利なライブラリを設計しようとしていると、この辺りを理解することはとても大事になります。 この辺りが一通り説明できたら、C++erレベル1は脱したと言ってもいいんじゃないかなと思います。

フィックスターズ広告

私は、実を言うと、ずいぶん前にstd::forwardを使えなくてあきらめたことがあったのですが、この発表の後、見返してみたら簡単に使えるようになってました。

フィックスターズは本物のC++プログラマーを募集しています
・・・というわけではないですが、高速化のお仕事が強いというだけあって、業務でもC++やC言語を使うことが多いですので、C++に強い関心がある方にもオススメです。

Work at Fixstars -フィックスターズで働いてみませんか?:

2016年2月19日金曜日

libSGM公開しました

NVIDIA Jetson プラットフォーム Meet-up #02 開催決定!

というわけで関係者の皆様、おめでとうございます。
僕も参加する予定です!(まだ登録してない
こ、今回はネタないかな…たぶん…

それはそうと、第一回で発表致しましたSGMですが、ライブラリとして公開しました。
というか、一週間前ぐらいからしていたのですが、なんかご連絡が遅くなりまして申し訳ありません。

名称はlibSGMとなっております。
公開はこちらにて行っております。
https://github.com/fixstars/libSGM

皆様是非使っていただいて、issueとかprとかお待ちしています。
あとこれ速くないぞというご指摘も歓迎致しますのでなにとぞなにとぞ…

さて、一週間前に公開したときは画像を入力としていたのですが
本日のコミットで、ステレオカメラのZEDをサポートし、ZEDで実行できるサンプルを同梱致しました。
Windows(+GPU環境)でも、JetsonTK1(Linux)でもZEDは動作すると思いますので、もしZEDお持ちの方は是非お試しいただけると嬉しいです。


ZEDでの動作風景が以下のような感じになっています。

この緑色の人は僕です。ものすごいアホっぽいポージングですね。
これだけだとちゃんと深度取れてるか怪しいので、ちゃんと深度が取れるようなポージングをしてみました。


さらにアホっぽくなりましたね。
ですが、突き出した手がちゃんと取れているように見えると思います。

ご質問等は遠慮なくgithubまでいただければ幸いです。
それでは皆様、ぜひ使ってくださいね!

2016年1月28日木曜日

NVIDIA Jetson プラットフォーム Meet-up #01で発表してきました

NVIDIA Jetson プラットフォーム Meet-up #01で発表してきました。

内容としては前回の記事を見ていただくほうがいいかもしれませんが、SGMチームで作っていたSGMについて、発表してきました。



Jetsonプラットフォーム Meet-upですが、非常に面白い会でした。
このような機会を設けていただいたNVIDIA様には、この場を借りて感謝を申し上げます。ありがとうございました。次もあるって期待しててもよろしいですか?

コード公開します、というお話をしたところ、結構興味を持っていただいて、非常に嬉しかったです。
まずはVisionWorksよりも速くすることを目的としたいと思います。

2016年1月18日月曜日

Realizing Self-Driving Cars with General-Purpose Processors 日本語版


Abstract: 自動運転システムというものが実用化されてきましたが、これはまだハードウェア実装に依存しています。しかしながら、これから数年のうちに、このようなコンピュートインテンシブな処理もソフトウェア実装で完璧にリアルタイム動作させることができるようになるでしょう。本稿では、私たちが組み込み用低消費電力GPU向けに実装・最適化したステレオ画像処理ベースのアルゴリズムについて述べます。私たちはこのシステムで、VGA画像処理で8FPSを達成しています。

Advanced Driver Assistance Systems(ADAS)と呼ばれる技術は、ここ数年、自動運転技術を実現するものとして、かなり注目を集めています。ダイムラー社が、2014年にメルセデス・ベンツEクラスとSクラスに搭載した6D-Vision Systemは、一つ大きな到達点となりました。
2013年時点で搭載されたものは、完全に自動運転というわけにはいきませんでしたが、Adaptive cruise control(日本語で言うと定速走行・車間距離制御装置)や、衝突回避システム、あるいはレーンキーピングシステム等を搭載しています。なかでも面白いのは、彼らがMAGIC BODY CONTROLと呼ぶ、Intelligentなサスペンションコントロールシステムです。
さて、もう少し最近で言えば、2015年5月に、ダイムラーグループのフレイトライナーが公開した、"Inspiration Truck"  があります。これは6D-Vision Systemをベースとしたもので、アメリカ合衆国はネバダ州の高速道路を走行することを法的に許可されたものです。
このシステムを構成しているのは演算量的に非常に"重たい"画像処理となっています。6D-Vision Systemのサイトによれば、このシステムのデータフローはFigure.1のようになります。



Figure 1: 6D-Vision Systemにおける画像処理の全体的なデータフロー
このシステムは、密なステレオカメラ処理――カメラ前方の各ピクセルごとの距離だとか――に依存するところが非常に大きいです。今まで様々な正確なアルゴリズムが考えだされてきましたが、2005年にSemi-Global Matching(SGM)と呼ばれるアルゴリズムが提案されるまで、"リアルタイム"なアプリケーション(VGA程度の大きさで30FPS以上、といったもの) を開発することはできませんでした。ダイムラーリサーチは2008年、市販のFPGAでこのアルゴリズムを使い、リアルタイムアプリケーションを開発しています。



Figure 2: ステレオ画像の例(Daimler Urban Segmentation Datasetより)

ステレオ処理のアルゴリズムというのは、それまで、計算量とロバスト性と正確性のバランスから、ブロックマッチングに代表されるアルゴリズムを使用していました。ほかには、動的計画法や、確率伝搬法、あるいはグラフカット問題などです。ですが、ロバスト性、正確性、計算効率のいずれをとっても、SGMと並ぶものは存在していません。SGMはブロックマッチングと比較すると計算量が多くはあるのですが、その正確性・ロバスト性は、計算コストの増加に見合う分の価値があります。実際、6D-VisionはSGMによるものです。



Figure 3: Figure2をOpenCVのブロックマッチングで視差を算出した画像(左)と、本稿で述べているSGMを使用して視差を算出した結果画像(右)
明るさが視差を表している。黒は不明領域


Figure 4: 左画像での座標(x, y)を中心とした矩形領域のブロックマッチングの例。この時、右画像でもっとも良いマッチング結果となるのは(x*, y)
従来のブロックマッチングでは、基準となるピクセルの周辺の領域(たとえば9x9の矩形とか) の画素値を比較することで、右画像と左画像の各画素を一致させようとします。この時の比較基準となる値をマッチングコストと呼称します。適切に画像を整えることにより、エピポーラ制約から、各ピクセルは同一のy座標を持つことが保証されるため、Figure 4に示す通り、x軸に沿って右の画像をスキャンするだけで済みます。ここで、点(x*, y)の周辺で、x軸におけるマッチングコストと、視差dが最小化する箇所を探索します。なぜなら、視差は深度情報に反比例するので、近くにある物体はより大きな視差を持つためです。このようにして得られた画像を、depth image(深度画像)と呼称します。


Figure 5: Stereo SGMが、座標(x, y)におけるスムージングコストを算出するために使用する8方向の図
しかしながら。Figure 3.を見てもらうとわかるように、まだ欠損や穴が見受けられます。それを軽減する方法の一つが、平滑化や、正規化といったものです。グラフカット問題や確率伝搬法などの従来の手法ではグローバルに平滑化しようとしますが、それではリアルタイムに良い結果を得ることはできません。一方で、"セミグローバルな"平滑化手法というのは革命的です。Figure 5.に示すように、グローバルな平滑化ではなく、各ピクセルの8方向程度に沿うことで、計算コストを最小限に抑えることができます。これはブロックマッチング同様、マッチングコストを最小化しているようですが、この場合では、道路やガードレール等の、十分な強度を得にくい、比較のしにくい表面においても行うことができ、"スムージングコスト"と呼称されます。キャリブレーションが正確に行われたステレオカメラで得られた三次元点群は、LIDAR(laser intensity direction and ranging、レーダー強度方向探知ならびに距離測定)で得られたものに匹敵します。


Figure 6: Segmentation SGMが座標(x,y)におけるスムージングコストを算出するための垂直方向の演算についての図
これで密な深度情報を得ることができました。次のステップとして、得られた情報を、道路とオブジェクトという二つの大きなクラスに分類します。道路が深度の大部分が水平に広がっている一方、オブジェクトは垂直に分布しています。しかし、Figure 3を見れば明らかなように、そのデプスマップは、カメラの傾きに依存して傾いています。分類を行う基本的なアイディアとして、まずメディアンフィルタを使用します。水平方向に5、または10ピクセルごとのメディアンフィルタを適用し、水平方向の画像幅の削減、ノイズの除去を行います。この水平方向のメディアンフィルタを適用している最中には、高さ方向への操作は行いません。その後、垂直方向にFigure 6のような、別のコスト最小化関数を適用します。総コスト関数は、マッチングコスト(測定深度と、オブジェクトや道路のモデル化された深度との差) だけでなく、モデル化された道路の事前情報等を持つ、スムージングコストとして提供しています。この状態で、道路かオブジェクトかの判定を行うには、すべてのピクセルにおいて、コスト関数の和を最小化します。コストの削減には、動的計画法がよく用いられてきましたが、SGMにもまた効果的に適用できることがわかりました。動的計画法は一次元の最適化問題に適しており、これはSGMが必要ないことを意味しています。しかし、SGMをステレオ、セグメンテーション(モデル分割) の双方に使用することで、ハードウェア・ソフトウェアを問わず、実装の手間を削減することができると信じています。


Figure 7: Segmentation SGMと、Figure2の画像を使用して得られた出力。
近いオブジェクトは赤く、遠いオブジェクトは緑色で表現されている。道路は透明。

現在、市販のFPGAによって、ハードウェアをカスタマイズするための労力は著しく下がっています。ここ数年から、あるいは今後数年かけて、ソフトウェアで書かれたコードをハードウェア化するソリューションが実用化されていくことでしょう。私たちの目的を証明するために、消費電力が10W以下のNVIDIA Tegra X1上でこのアルゴリズムを実装し、最適化を行いました。その処理速度がTable 1.になります。NVIDIAのロードマップに依れば、2020年前後で、SGMの二つのアルゴリズム(Stereo SGM, Segmentation SGM) をVGAサイズでリアルタイムに動作させることができるようになると思われます。
 組み込み分野に対して(フィックスターズ含む)、新規の参入者ができることは、高い計算コストを持つアルゴリズムをソフトウェアで完全に実装し、それをADASを開発している人々に提供していくことであると考えています。

Width x Height x Max DisparityAverage processing time (ms)
640x480x641024x440x128
Stereo SGM80282
Segmentation SGM3942
Table 1: Tegra X1上での、Stereo SGM、Segmentation SGMの
現在のパフォーマンス
さて、デモムービーを用意しました。この動画の左右画像は、ダイムラー社が提供しているDaimler Urban Segmentation Dataset.を利用しています。わかりやすくするために、セグメント化されたオブジェクトには輪郭処理を、距離が接近したオブジェクトには警告マークを付与して表示しています。また、このデモではセグメンテーションの精度よりも速度を優先しています。もともとのSegmentation SGMのアルゴリズムでは、ノイズは少ないのですが、並列化が非常に困難です。そのため、まだ改善の余地があると考えています。現状では、オブジェクトのトラッキングを行わないため、移動している物体が自分と衝突するコースにいるかどうかという点を知ることはできません。

オブジェクトを追跡する前に、まず、オプティカルフローの密な推定値を取得する必要があります。伝統的には、ロバスト性と正確性を両立させるために、二つのアルゴリズムを使用します。特徴点と特徴量に基づいたグローバルでスパースな構造(CENSUS, SIFT, etc) と、局所改善の二つです。これらの代わりに、私たちはflow SGMというアルゴリズムを使用して、この問題に取り組み始めています。flow SGMはSegmentation SGMと同様に、Stereo SGMの派生として知られています。この手法によって、従来と同等の結果を得つつ、Stereo SGM, Segmentation SGMからのコードを再利用することで、システム全体の複雑さを軽減できることを願っています。現在、予備実験の結果は良いものが得られているため、今後、またブログにて調査結果を公開する予定です。

- The "SGM Team": Samuel Audet, Yoriyuki Kitta, Yuta Noto, Ryo Sakamoto, Akihiro Takagi

About us

About us
このブログは、株式会社フィックスターズのエンジニアが、あらゆるテーマについて自由に書いているブログです。

Popular Posts

Powered by Blogger.
Subscribe Via Email

Subscribe to our newsletter to get the latest updates to your inbox. ;-)

Your email address is safe with us!