未分類

CUDAで二次元配列の転置は高効率にできないと言ったな！あれは嘘だ！

by belltail
1
寒くなってくるとお魚で鍋とかもしたくなってきますね．
海辺のお魚屋さんでは人々が行列を成してるんでしょうか．

行列といえばGPGPUが活躍する最も重要な分野の一つと言えますね．

んで，そんな行列でよく使われるタスクのひとつに「転置」があります．

\bold{A}= \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \\ \end{pmatrix}, \bold{A}^{\mathrm{T}}= \begin{pmatrix} 1 & 4 & 7 \\ 2 & 5 & 8 \\ 3 & 6 & 9 \\ \end{pmatrix}
\bold{B}= \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ \end{pmatrix}, \bold{B}^{\mathrm{T}}= \begin{pmatrix} 1 & 4 \\ 2 & 5 \\ 3 & 6 \\ \end{pmatrix}

画像屋さんなら，線形なフィルタをかけるときには縦方向にフィルタして横方向にフィルタすればいいのですが，実装的には横方向にフィルタして画像全体を転置してまた横方向にフィルタして転置して戻す方法のほうが効率的なことが多いです．特にGPGPUなどでは．

ということでよく使う転置というタスクですが，CUDAでやるにはちょっと工夫しないとメモリアクセス効率が下がります．今日のお題はそれで．

  1. CUDAでいちばん単純な転置

    1. データの配置
    2. ということは
    3. 単純にやるとメモリアクセス効率が悪い
  2. 改善策

    1. 改良したカーネル
    2. shared memoryのbank conflictには注意
  3. 実行速度まとめ
  4. 全体コード

    1. main.cpp
    2. kernel.cu
    3. Makefile
  5. まとめ

CUDAでいちばん単純な転置

データの配置

行オーダーで記録された画像があるとします．つまり，
I(x,y)=\left[ \begin{array}{ccccc} I(0,0) & I(1,0) & I(2,0) & & I(w - 1, 0) \\ I(0,1) & I(1,1) & I(2,1) & \cdots & I(w - 1, 0) \\ I(0,2) & I(1,2) & I(2,2) & & I(w - 1, 0) \\ & \vdots & & \ddots & \vdots \\ I(0,h-1) & I(1,h-1) & I(2,h-1) & \ldots & I(w - 1, h-1) \\ \end{array} \right]
な画像もしくは行列があるとして，これが0行目をだーっと並べた後そのまま1行目をだーっと並べて…という形でw \times hの一次元配列に収まる形です．

CUDAの特性に併せて128bytes境界に合わせるためのパディングをする，などは考慮しません（今回の場合はそれをやっても効果はないと思います）．

ということは

座標x,yの，配列中における位置はそのまんまy \times w+xで，転置先はx \times h + yになりますね．それをそのまま実装したカーネルが以下のようになります．

template<class T>
__global__ void transpose(T* d_dst, const T* d_src, int w, int h)
{
    const int sx = blockDim.x * blockIdx.x + threadIdx.x;
    const int sy = blockDim.y * blockIdx.y + threadIdx.y;
    if(w <= sx || h <= sy)
         return;
    d_dst[sx * h + sy] = d_src[sy * w + sx];
}
単純にやるとメモリアクセス効率が悪い

たとえばフルHDのグレースケール画像を転置してみると，カーネルの実行に0.257msかかりました．

この方法だと，読み込みのアクセスはCoalescedされていますが，書き込みのアクセスがCoalescedされていません．
これだけで意味がわかる方は次の説明は飛ばしてしまってOKです．

まずCUDA一般論的な話として，Coalesced Accessとはなにか，詳細は他に譲るとして簡単に説明．

  • CUDAではWarp=32threadsを単位としてSIMT実行される(Fermi系以降？確か)
  • メモリアクセスは128bytes境界でまとめてアクセスされる
  • なので，warp内の各スレッドがメモリ中の連続領域に同時アクセスすると，メモリアクセスは1回ですむ

    • 厳密には，warp内の各スレッドが同時にアクセスする先が全て128bytes境界の中であれば1回にまとめられる
    • 128bytes境界とは，簡単に言うとアドレスが128の倍数で始まる128バイトの領域区間
  • warp内の各スレッドがメモリ中のばらばらの領域に同時アクセスすると，それぞれのアクセス先でそこを含む128bytes境界単位でアクセスされる
  • なのでメモリアクセスパターンが悪いと，実効メモリ速度が理論上32分の1にまで低下する

    • GTX 580やGTX 680の場合は，6GB/sが上限になるわけですね…

これを上のカーネルに当てはめると，

実際にVisual Profilerで測定してみると，

メモリアクセス効率がかなり悪いんですね（特にGlobal Store Efficiency）

改善策

要は，各warpで横一列に取ってきたデータを，縦一列で書きこまないといけないことが問題なのです．
横一列で読んで横一列で書くことができれば，連続領域への書き込みになってアクセス効率が上がります．

イメージとしてはこんなかんじ．

もちろんそのままではできなくて，

みたいにすればいいんじゃね？というのがアイデア．はい，図じゃわかりません＞＜

具体的には，shared memoryを使うことで，読んだ画素と書く画素を変えることで実現できます．

改良したカーネル
template<class T, int BLOCK_DIM>
__global__ void transpose_shared(T* d_dst, const T* d_src, int w, int h)
{
    __shared__ T buf[BLOCK_DIM][BLOCK_DIM + 1];   //bank conflictを避けるための+1
    const int sx = blockDim.x * blockIdx.x + threadIdx.x;
    const int sy = blockDim.y * blockIdx.y + threadIdx.y;
    if(w <= sx || h <= sy)
        return;

    buf[threadIdx.x][threadIdx.y] = d_src[sy * w + sx];
    __syncthreads();

    d_dst[(blockIdx.x * blockDim.x + threadIdx.y) * h + (blockIdx.y * blockDim.y + threadIdx.x)] = buf[threadIdx.y][threadIdx.x];
}

注目箇所は，d_dstに書きこむ時のアドレッシングです．
どーなってんのかよくわからない＞＜という方は上の図ともにらめっこしながら・・・

ひとまずこれで，Visual ProfilerでGlobal Memoryの効率を見ると

Store Efficiencyが100%になりました(・∀・)
Loadが低いのはなんでですかね？（を

shared memoryのbank conflictには注意

shared memoryは，4バイトごとに32のバンク（Fermi以降）に分かれてています．なので，8×8や16×16の正方形のshared memoryを確保してしまうと激しくbank conflictします．
そのため1要素だけずらすことでこれを回避しています．

何も考慮しない場合，shared memoryへのアクセスに占めるconflictの割合を表すshared memory replay overheadが高いです．

これを，こんなふうに減らすことができます．

実行速度まとめ

グレースケールのフルHDを処理したときの，カーネル実行時間

デバイス カーネル 実行時間 レガシー比速度
GTX 670 レガシー 0.257ms +0%
GTX 670 shared (low conflict) 0.216ms +19%
GTX 670 shared (high conflict) 0.236ms +9%
GTX 580 レガシー 0.228ms +0%
GTX 580 shared (low conflict) 0.147ms +55%
GTX 580 shared (high conflict) 0.183 +25%

というかんじになりました．

GTX 670では，そもそも速度差があんまりなくて，conflictの影響も比較的少ない．
GTX 580では，かなり効果は顕著で，conflictもけっこう効く．と言ったところでしょうか．

GTX 670ってFLOPSではGTX 580とは比べ物にならないほど速いし，メモリ速度も192GB/sで同じなのに，実アプリでは逆転してこの差．悲しい・・・．

全体コード

CUDA 4.1以降とOpenCV 2.3以降を想定してます．
時間計測の関係でカーネルの終了を同期的に待ったり，関数を細かく分けたりとかごちゃごちゃしたことをしてます．

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

void initialize(int _w, int _h);
void send(unsigned *h_src);
void transpose();
void restore(unsigned *h_dst);

cv::Mat myTranspose(const cv::Mat& img)
{
    const int w = img.cols;
    const int h = img.rows;

    std::vector<unsigned> buf(w * h);
    for(int y = 0; y < h; ++y)
        for(int x = 0; x < w; ++x)
            buf[y * w + x] = static_cast<unsigned>(img.at<unsigned char>(y, x));

    initialize(w, h);
    send(&buf[0]);
    transpose();
    restore(&buf[0]);

    cv::Mat img2(cv::Size(h, w), CV_8UC1);
    for(int y = 0; y < w; ++y)
        for(int x = 0; x < h; ++x)
            img2.at<unsigned char>(y, x) = static_cast<unsigned char>(buf[y * h + x]);
    return img2;
}

int main(int argc, char *argv[])
{
    //モノクロでロード
    cv::Mat img = cv::imread(argv[1], 0);
    cv::Mat img2 = myTranspose(img);
    cv::imwrite("result.png", img2);
}

kernel.cu
dim3 calcBlock(dim3 thread, int w, int h)
{
    return dim3(
        static_cast<int>(std::ceil(1.0 * w / thread.x)),
        static_cast<int>(std::ceil(1.0 * h / thread.y)));
}

template<class T>
__global__ void transpose(T* d_dst, const T* d_src, int w, int h)
{
    const int sx = blockDim.x * blockIdx.x + threadIdx.x;
    const int sy = blockDim.y * blockIdx.y + threadIdx.y;
    if(w <= sx || h <= sy)
        return;
    d_dst[sx * h + sy] = d_src[sy * w + sx];
}

template<class T, int BLOCK_DIM>
__global__ void transpose_shared(T* d_dst, const T* d_src, int w, int h)
{
    __shared__ T buf[BLOCK_DIM][BLOCK_DIM + 1];
    const int sx = blockDim.x * blockIdx.x + threadIdx.x;
    const int sy = blockDim.y * blockIdx.y + threadIdx.y;
    if(w <= sx || h <= sy)
        return;

    buf[threadIdx.x][threadIdx.y] = d_src[sy * w + sx];
    __syncthreads();

    d_dst[(blockIdx.x * blockDim.x + threadIdx.y) * h + (blockIdx.y * blockDim.y + threadIdx.x)] = buf[threadIdx.y][threadIdx.x];
}

int w, h;
unsigned* d_src;
unsigned* d_dst;

void initialize(int _w, int _h)
{
    w = _w;
    h = _h;
    //領域確保
    cudaMalloc((void**)&d_src, sizeof(unsigned) * w * h);
    cudaMalloc((void**)&d_dst, sizeof(unsigned) * h * w);
    cudaThreadSynchronize();
}

void send(unsigned *h_src)
{
    //画像の転送
    cudaMemcpy(d_src, h_src, sizeof(unsigned) * w * h, cudaMemcpyHostToDevice);
    cudaThreadSynchronize();
}

//転置処理
void transpose()
{
    dim3 thread(16, 16);
    //ここでカーネルを切り替え
    //transpose<<<calcBlock(thread, w, h), thread>>>(d_dst, d_src, w, h);
    transpose_shared<unsigned, 16><<<calcBlock(thread, w, h), thread>>>(d_dst, d_src, w, h);
    cudaThreadSynchronize();
}

void restore(unsigned *h_dst)
{
    //結果を書き戻す
    cudaMemcpy(h_dst, d_dst, sizeof(unsigned) * w * h, cudaMemcpyDeviceToHost);
    cudaThreadSynchronize();
}
Makefile

ついでにミニマムなMakefileだけかいときました

TARGET		= prog.out
NVCC		= nvcc
SOURCE		= main.cpp
CUDA_SOURCE	= kernel.cu
CXXFLAGS	= -I$(HOME)/include -std=c++0x -O2
CUDAFLAGS	= 
LDFLAGS		= -L$(HOME)/lib -lopencv_core -lopencv_highgui -lopencv_imgproc

$(TARGET): $(SOURCE) $(CUDA_SOURCE)
	$(CXX) -c $(SOURCE) $(CXXFLAGS)
	$(NVCC) -c $(CUDA_SOURCE) $(CUDAFLAGS)
	$(NVCC) *.o $(LDFLAGS) -o $(TARGET)

clean:
	$(RM) $(TARGET) *.o

まとめ

shared memoryのいちばん大きな用途は「制御されたキャッシュ」としてのものだと思いますが，「ブロック内スレッド間でのデータ共有」ができることを活用してこのように「アクセスパターン最適化」のためにも活かすことができます．
うまく使って速いコード書きましょう．

ついでに，shared memoryでのブロック内データ共有もいいですが，Keplerから使える「shuffle」を使うと一部代替できるかもしれません．
これが使える局面であれば，conflictのことを気にする必要がなくオーバヘッドも少ないのでさらに高性能になります．

One thought on “CUDAで二次元配列の転置は高効率にできないと言ったな！あれは嘘だ！”

  1. ピンバック: Words for GPGPU | STUDY DIARY

コメントを残す

メールアドレスが公開されることはありません。 ※ が付いている欄は必須項目です