10

aとbの2つの配列があり、「最小畳み込み」を計算して結果cを生成したいと思います。単純な擬似コードは次のようになります。

for i = 0 to size(a)+size(b)
    c[i] = inf
    for j = 0 to size(a)
        if (i - j >= 0) and (i - j < size(b))
            c[i] = min(c[i], a[j] + b[i-j])

（編集：ループを1ではなく0から開始するように変更）

minが代わりに合計である場合、高速フーリエ変換（FFT）を使用できますが、minの場合、そのようなアナログはありません。代わりに、GPU（CUDA）を使用して、この単純なアルゴリズムをできるだけ高速にしたいと思います。これを実行する既存のコード（またはFFTなしで合計ケースを実装するコードを見つけて、目的に合わせて調整できるようにする）を見つけてうれしいですが、これまでのところ、良い結果は得られていません。私のユースケースには、サイズが1,000〜100,000のaとbが含まれます。

質問：

  • これを効率的に行うためのコードはすでに存在しますか？

  • これを自分で実装する場合、構造的に、効率を最大化するためにCUDAカーネルはどのように見えるべきですか？各c[i]が別々のスレッドによって計算される単純なソリューションを試しましたが、これは最善の方法ではないようです。スレッドブロック構造とメモリアクセスパターンを設定する方法に関するヒントはありますか？

4

3 に答える 3

5

a大規模な場合に役立つ代替手段は、出力エントリごとにブロックbを使用することです。ブロックを使用すると、メモリの結合が可能になります。これは、メモリ帯域幅が制限された操作で重要になります。また、かなり効率的な共有メモリの削減を使用して、スレッドごとの部分的な結果をブロックごとの最終的な結果に結合することができます。おそらく最良の戦略は、MP ごとに同時に実行されるのと同じ数のブロックを起動し、各ブロックが複数の出力ポイントを発行するようにすることです。これにより、総命令数が比較的少ない多くのブロックの起動とリタイアに関連するスケジューリング オーバーヘッドの一部が解消されます。c

これがどのように行われるかの例:

#include <math.h>

template<int bsz>
__global__ __launch_bounds__(512)
void minconv(const float *a, int sizea, const float *b, int sizeb, float *c)
{
    __shared__ volatile float buff[bsz];
    for(int i = blockIdx.x; i<(sizea + sizeb); i+=(gridDim.x*blockDim.x)) {
        float cval = INFINITY;
        for(int j=threadIdx.x; j<sizea; j+= blockDim.x) {
            int t = i - j;
            if ((t>=0) && (t<sizeb))
                cval = min(cval, a[j] + b[t]);
        }
        buff[threadIdx.x] = cval; __syncthreads();
        if (bsz > 256) {
            if (threadIdx.x < 256) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+256]);
            __syncthreads();
        }
        if (bsz > 128) {
            if (threadIdx.x < 128) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+128]); 
            __syncthreads();
        }
        if (bsz > 64) {
            if (threadIdx.x < 64) 
                buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+64]);
            __syncthreads();
        }
        if (threadIdx.x < 32) {
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+32]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+16]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+8]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+4]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+2]);
            buff[threadIdx.x] = min(buff[threadIdx.x], buff[threadIdx.x+1]);
            if (threadIdx.x == 0) c[i] = buff[0];
        }
    }
}

// Instances for all valid block sizes.
template __global__ void minconv<64>(const float *, int, const float *, int, float *);
template __global__ void minconv<128>(const float *, int, const float *, int, float *);
template __global__ void minconv<256>(const float *, int, const float *, int, float *);
template __global__ void minconv<512>(const float *, int, const float *, int, float *);

[免責事項: テストもベンチマークもされていません。自己責任で使用してください]

これは単精度浮動小数点ですが、同じ考え方が倍精度浮動小数点でも機能するはずです。INFINITY整数の場合、C99マクロをINT_MAXまたはのようなものに置き換える必要がありますLONG_MAXが、原則はそれ以外は同じです。

于 2012-11-06T22:26:21.123 に答える
5

より高速なバージョン:

__global__ void convAgB(double *a, double *b, double *c, int sa, int sb)
{
    int i = (threadIdx.x + blockIdx.x * blockDim.x);
    int idT = threadIdx.x;
    int out,j;

    __shared__ double c_local [512];

    c_local[idT] = c[i];

    out = (i > sa) ? sa : i + 1;
    j   = (i > sb) ? i - sb + 1 : 1;

    for(; j < out; j++)
    {    
       if(c_local[idT] > a[j] + b[i-j])
          c_local[idT] = a[j] + b[i-j]; 
    }   

    c[i] = c_local[idT];
} 

**Benckmark:**
Size A Size B Size C Time (s)
1000   1000   2000   0.0008
10k    10k    20k    0.0051
100k   100k   200k   0.3436
1M     1M     1M     43,327

古いバージョン、1000 から 100000 の間のサイズの場合、この素朴なバージョンでテストしました:

__global__ void convAgB(double *a, double *b, double *c, int sa, int sb)
{
    int size = sa+sb;

    int idT = (threadIdx.x + blockIdx.x * blockDim.x);
    int out,j;


    for(int i = idT; i < size; i += blockDim.x * gridDim.x)
    {
        if(i > sa) out = sa;
        else out = i + 1;

        if(i > sb) j = i - sb + 1;
        else j = 1;


        for(; j < out; j++)
        {
                if(c[i] > a[j] + b[i-j])
                    c[i] = a[j] + b[i-j];
        }
    }
}

配列にいくつかのランダムな倍精度数と 999999 (テスト用) を入力aしbましcた。c関数を使用して（変更なしで）（CPU内の）配列を検証しました。

また、内側のループの内側から条件を削除したので、一度だけテストします。

100% 確信があるわけではありませんが、次の変更は理にかなっていると思います。を持っていたのでi - j >= 0、これは と同じで、これは、このブロック 'X' に入らないとi >= jすぐに(j++ から) ということを意味します。j > i

if(c[i] > a[j] + b[i-j])
   c[i] = a[j] + b[i-j];

したがって、変数outでループ条件 ifを計算しましたi > sa。これは、ループが終了するときにループが終了することをj == sa意味i < saしi + 1ますi >= j。

他の条件は、i - j < size(b)ブロック 'X' の実行を開始することを意味します。i > size(b) + 1jj

if(i > sb) j = i - sb + 1;
else j = 1;

このバージョンを実際のデータ配列でテストできるかどうかを確認し、フィードバックをお寄せください。また、改善点は大歓迎です。

編集:新しい最適化を実装できますが、これは大きな違いはありません。

if(c[i] > a[j] + b[i-j])
    c[i] = a[j] + b[i-j];

次の方法でifを排除できます。

double add;
...

 for(; j < out; j++)
 {
   add = a[j] + b[i-j];
   c[i] = (c[i] < add) * c[i] + (add <= c[i]) * add;
 }

持つ：

if(a > b) c = b; 
else c = a;

c = (a < b) * a + (b <= a) * b と同じです。

a > b の場合、c = 0 * a + 1 * b; => c = b; a <= b の場合、c = 1*a + 0 *b; => c = a;

**Benckmark:**
Size A Size B Size C Time (s)
1000   1000   2000   0.0013
10k    10k    20k    0.0051
100k   100k   200k   0.4436
1M     1M     1M     47,327

CPU から GPU へのコピー、カーネルの実行、GPU から CPU へのコピーの時間を測定しています。

GPU Specifications   
Device                       Tesla C2050
CUDA Capability Major/Minor  2.0
Global Memory                2687 MB
Cores                        448 CUDA Cores
Warp size                    32
于 2012-11-06T00:48:23.393 に答える
2

私はあなたのアルゴリズムを使用しました。お役に立てると思います。

const int Length=1000;

__global__ void OneD(float *Ad,float *Bd,float *Cd){
    int i=blockIdx.x;
    int j=threadIdx.x;
    Cd[i]=99999.99;
    for(int k=0;k<Length/500;k++){
        while(((i-j)>=0)&&(i-j<Length)&&Cd[i+k*Length]>Ad[j+k*Length]+Bd[i-j]){
            Cd[i+k*Length]=Ad[j+k*Length]+Bd[i-j];
    }}}

ブロックごと500にスレッドを取得しました。そして、グリッドごとのブロック。私のデバイスのブロックあたりのスレッド数は に制限されているため、スレッドを使用しました。すべての配列のサイズを(=1000) としました。500512500Length

働く：

  1. iブロック インデックスをj 格納し、スレッド インデックスを格納します。

  2. スレッドの for数が配列のサイズよりも少ないため、ループが使用されます。

  3. while ループは繰り返しに使用されますCd[n]。

  4. 多くのブロックとスレッドを使用したため、共有メモリを使用していません。したがって、各ブロックに必要な共有メモリの量は少なくなります。

PS:デバイスがより多くのスレッドとブロックをサポートしている場合はk<Length/500、k<Length/(supported number of threads)

于 2012-11-07T18:58:43.933 に答える