So-net無料ブログ作成
検索選択

CUDAで離散ウェーブレット変換 [プログラミング]

CUDAで離散ウェーブレット変換を行って見ましょう。

離散ウェーブレット変換には、一番簡単な、Haar wavelet を使います。

とりあえず、作ってみました。
アルゴリズムは、Wikipedia の 離散ウェーブレット変換 のページ にあるものを参考にしました。
また、以前作ったメモリクラスを使用しています。

#include  <cutil_inline.h>
#include  <iostream>
#include  "cuda_mem.h"

__global__ void wavelet(const float* input, float* output, int N)
{
    // 前提:
    // N = 2^m = blockDim.x * 2
    // sizeof (shared) = sizeof(float) * N * 3

    extern __shared__ float shared[];

    float* in = shared;        // 入力データ
    float* buff = shared + N;  // 中間データ
    float* out = shared + N*2; // 計算結果

    // 入力データをシェアードメモリへコピー
    in[threadIdx.x] = input[threadIdx.x];
    in[threadIdx.x + blockDim.x] = input[threadIdx.x + blockDim.x];
    __syncthreads();

    // 計算
    for (int n = N / 2; n >= 1; n /= 2) {
        if (threadIdx.x < n) {
            buff[threadIdx.x] = in[threadIdx.x * 2] + in[threadIdx.x * 2 + 1];
            out[n + threadIdx.x] = in[threadIdx.x * 2] - in[threadIdx.x * 2 + 1];
        }

        // 入力データと中間データを入れ替え
        float* temp = in;
        in = buff;
        buff = temp;
        __syncthreads();
    }

    if (threadIdx.x == 0)
        out[0] = in[0];

    // 計算結果をグローバルメモリへコピー
    output[threadIdx.x] = out[threadIdx.x];
    output[threadIdx.x + blockDim.x] = out[threadIdx.x + blockDim.x];
    __syncthreads();
}

int main()
{
    const int N = 1024;        // データ数。1024 以下の2のべき乗

    CudaMem<float> data(N);
    CudaMem<float> result(N);

    for (int i = 0; i < N; i++)  // 入力データ取得
        std::cin >> data[i];

    data.to_gpu();
    wavelet<<<1, N/2, sizeof(float)*N*3>>>(data.get_gpu(), result.get_gpu(), N);

    result.from_gpu();

    for (int i = 0; i < N; i++)
        std::cout << result[i] << "\n";
    
    return 0;
}

複数のブロックを使うと色々と面倒そうだったので、1つのブロックで処理できるようデータ数の上限は 1024 です。
(ブロック内のスレッド数上限が 512 だし、シェアードメモリのサイズ上限が16kbyte)


ところで、上のプログラムでは、

            buff[threadIdx.x] = in[threadIdx.x * 2] + in[threadIdx.x * 2 + 1];
            out[n + threadIdx.x] = in[threadIdx.x * 2] - in[threadIdx.x * 2 + 1];

この部分で、シェアードメモリのバンク コンフリクトが発生しています。
(バンク コンフリクトが発生するとシェアードメモリへのアクセスに若干のロスが発生する)

バンク コンフリクトを解消するには、
シェアードメモリのアクセスが、shared_mem[threadIdx.x * A + B] という形式の場合、A が奇数になるようにします。

そこで、データの位置を並べ替え、こんな↓風に処理できるようにします。

            buff[threadIdx.x] = in[threadIdx.x] + in[n + threadIdx.x];
            out[n + threadIdx.x] = in[threadIdx.x] - in[n + threadIdx.x];

計算結果が変わらないようにするには、次の表のようにデータを並べ替えます。(データ数 = 8 の場合)

01234567
data[0]data[4]data[2]data[6]data[1]data[5]data[3]data[7]

一見、どう並べ変えているのか分かりづらいのですが、配列の添え字を 2進数で表現すると気付きます。
添え字を2進数で表してビットの並びを逆順にします。

data[13] の場合、13 は2進数で 1101。
これを逆順にして 1011 (データ数が16個の場合)。
よって、data[13]は、配列の11番目に移動します。

ちなみにデータ数が32個の場合、13 は 01101 と表して、この逆順は 10110 になります。

あと、離散ウェーブレット変換した結果も順番がばらばらになっているので正しい順番に戻す必要がありますが、これは、計算結果配列の添え字を2進数で表し、1になっている一番上位の桁を残して、それより下位の桁を逆順にします。
1011 は、一番左の 1 はそのままで、残りの 011 を逆順にして 1110 にします。

そんなわけで書き換えてみました。

#include  <cutil_inline.h>
#include  <iostream>
#include  "cuda_mem.h"

__global__ void wavelet(const float* input, float* output, int N)
{
    // 前提:
    // N = 2^m = blockDim.x * 2
    // sizeof (shared) = sizeof(float) * N * 3

    extern __shared__ float shared[];

    float* in = shared;        // 入力データ
    float* buff = shared + N;  // 中間データ
    float* out = shared + N*2; // 計算結果

    // 入力データをシェアードメモリへコピー (後の処理のため順番入れ替え)
    in[reverse_bit(threadIdx.x, M)] = input[threadIdx.x];
    in[reverse_bit(threadIdx.x + blockDim.x, M)] = input[threadIdx.x + blockDim.x];
    __syncthreads();

    // 計算
    for (int n = N / 2, m = N - 1; n >= 1; n /= 2, m--) {
        if (threadIdx.x < n) {
            buff[threadIdx.x] = in[threadIdx.x] + in[n + threadIdx.x];
            out[n + reverse_bit(threadIdx.x, m)] = in[threadIdx.x] - in[n + threadIdx.x];
        }

        // 入力データと中間データを入れ替え
        float* temp = in;
        in = buff;
        buff = temp;
        __syncthreads();
    }

    if (threadIdx.x == 0)
        out[0] = in[0];

    // 計算結果をグローバルメモリへコピー
    output[threadIdx.x] = out[threadIdx.x];
    output[threadIdx.x + blockDim.x] = out[threadIdx.x + blockDim.x];
    __syncthreads();
}

// value の下位 bitビットの並びを逆転させる
__device__ int reverse_bit(int value, int bit)
{
    略
}

int main()
{
    const int M = 10;
    const int N = 1 << M;    // データ数。1024 以下の2のべき乗

    CudaMem<float> data(N);
    CudaMem<float> result(N);

    for (int i = 0; i < N; i++)  // 入力データ取得
        std::cin >> data[i];

    data.to_gpu();
    wavelet<<<1, N/2, sizeof(float)*N*3>>>(data.get_gpu(), result.get_gpu(), N, M);

    result.from_gpu();

    for (int i = 0; i < N; i++)
        std::cout << result[i] << "\n";
    
    return 0;
}

ここで気付きました。

バンク コンフリクトによるロスより、ビットの逆転処理×2の方が時間がかかるんじゃないかと。

そして、分かりづらいのですが、データの並べ替えをしているところでバンク コンフリクトが発生しているんじゃないかと。
CUDA SDK に含まれている {CUDA SDKのインストール先}\C\common\bank_checker.h が バンク コンフリクトが発生しているかをチェックするクラスなので、これを使って確認するのも良いかも。使い方が分からないけど。


タグ:Cuda
nice!(0)  コメント(0)  トラックバック(0) 

nice! 0

コメント 0

コメントを書く

お名前:
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

トラックバック 0

この記事のトラックバックURL:
※言及リンクのないトラックバックは受信されません。

この広告は前回の更新から一定期間経過したブログに表示されています。更新すると自動で解除されます。