17

私は C++ と CUDA/C を使用しており、特定の問題のコードを書きたいと思っています。非常にトリッキーなリダクションの問題に遭遇しました。

並列プログラミングでの私の経験は無視できるものではありませんが、非常に限られており、この問題の特異性を完全に予測することはできません。私が直面している問題を処理するための便利な、または「簡単な」方法があるとは思えませんが、おそらく私は間違っています。この問題または同様の問題をカバーするリソース (記事、書籍、Web リンクなど) またはキーワードがある場合は、お知らせください。

あまりにも多くのコードを投稿するのではなく、ケース全体を可能な限り一般化し、抽象化するようにしました。

レイアウト ...

N 個の初期要素と N 個の結果要素のシステムがあります。(たとえば、N=8 を使用しますが、N は 3 より大きい任意の整数値にすることができます。)

static size_t const N = 8;
double init_values[N], result[N];

自己干渉なしで、init-values のほぼすべての (すべてではありませんが) 一意の順列を計算する必要があります。

これは、計算f(init_values[0],init_values[1])、f(init_values[0],init_values[2])、 ...、f(init_values[0],init_values[N-1])、f(init_values[1],init_values[2])、 ... 、 ... などを意味f(init_values[1],init_values[N-1])します。

これは、実際には、次の図に示す形状の仮想三角行列です。

 P     0    1    2    3    4    5    6    7
   |---------------------------------------
  0|   x 
   |
  1|   0    x
   |
  2|   1    2    x 
   |
  3|   3    4    5    x
   |
  4|   6    7    8    9    x
   |
  5|  10   11   12   13   14    x
   |
  6|  15   16   17   18   19   20    x
   |
  7|  21   22   23   24   25   26   27    x

各要素は、 のそれぞれの列要素と行要素の関数ですinit_values。

P[i] (= P[row(i)][col(i]) = f(init_values[col(i)], init_values[row(i)])

すなわち

P[11] (= P[5][1]) = f(init_values[1], init_values[5])

例 を使用すると、可能な一意の組み合わせがあります(N*N-N)/2 = 28(注: P[1][5]==P[5][1]、したがって、下 (または上) 三角行列しかありません) N = 8。

基本的な問題

結果の配列は、行要素の合計からそれぞれの列要素の合計を差し引いたものとして P から計算されます。たとえば、位置 3 の結果は、行 3 の合計から列 3 の合計を引いて計算されます。

result[3] = (P[3]+P[4]+P[5]) - (P[9]+P[13]+P[18]+P[24])
result[3] = sum_elements_row(3) - sum_elements_column(3)

N=4の絵で説明してみました。

必要な三角縮小スキーム

結果として、次のことが当てはまります。

  • N-1操作 (潜在的な同時書き込み) はそれぞれに対して実行されますresult[i]
  • result[i]減算と加算N-(i+1)からの書き込みがありますi
  • それぞれからの発信P[i][j]には、への減算とへの加算がありr[j]ますr[i]

ここで主な問題が発生します。

  • 1 つのスレッドを使用して各 P を計算し、結果を直接更新すると、複数のカーネルが同じ結果の場所 (それぞれ N-1 スレッド) に書き込もうとすることになります。
  • 一方、後続の削減ステップのために行列 P 全体を保存することは、メモリ消費の点で非常に高価であり、したがって、非常に大規模なシステムでは不可能です。

スレッドブロックごとに一意の共有結果ベクトルを持つという考えも不可能です。(N of 50k は 25 億の P 要素を作成するため、[ブロックあたりの最大数が 1024 スレッドであると仮定して] 各ブロックに 50k の倍精度要素を持つ独自の結果配列がある場合、最小数の 240 万ブロックが 900GiB を超えるメモリを消費します。)

より静的な動作の削減を処理できると思いますが、この問題は潜在的な同時メモリ書き込みアクセスに関してかなり動的です。（または、「基本的な」タイプの削減で処理できますか？）

いくつかの複雑さを追加します...

残念ながら、初期値に依存しない (任意のユーザー) 入力によっては、P の一部の要素をスキップする必要があります。順列 P[6]、P[14]、および P[18] をスキップする必要があるとします。したがって、計算する必要がある 24 の組み合わせが残っています。

どの値をスキップする必要があるかをカーネルに伝える方法は? N が非常に大きい場合 (数万の要素など)、それぞれに顕著な欠点があります。

1.すべての組み合わせを保存...

...それぞれの行と列の indexを使用して、 a で計算し、このベクトルで操作するstruct combo { size_t row,col; };必要があります。vector<combo>(現在の実装で使用)

std::vector<combo> elements;
// somehow fill
size_t const M = elements.size();
for (size_t i=0; i<M; ++i)
{
    // do the necessary computations using elements[i].row and elements[i].col  
}

このソリューションは、「数」（数万の要素でさえあるかもしれませんが、合計で数十億とは対照的ではありません）であるため、大量のメモリを消費しますが、回避します

  • 指数計算
  • 削除された要素の検索

これは、2 番目のアプローチの欠点です。

2. P のすべての要素を操作し、削除された要素を見つける

P の各要素を操作し、ネストされたループ (cuda ではうまく再現できませんでした) を回避したい場合は、次のようにする必要があります。

size_t M = (N*N-N)/2;
for (size_t i=0; i<M; ++i)
{
   // calculate row indices from `i`
   double tmp = sqrt(8.0*double(i+1))/2.0 + 0.5;
   double row_d = floor(tmp);
   size_t current_row = size_t(row_d);
   size_t current_col = size_t(floor(row_d*(ict-row_d)-0.5));
   // check whether the current combo of row and col is not to be removed
   if (!removes[current_row].exists(current_col))
   {
     // do the necessary computations using current_row and current_col
   }
}

このベクトルremovesは、最初の例のベクトルとは対照的に非常に小さいですが、elementsを取得するための追加の計算と if 分岐は非常に非効率的です。(まだ数十億の評価について話していることを思い出してください。)current_rowcurrent_col

3. P のすべての要素を操作し、後で要素を削除する

私が持っていた別のアイデアは、すべての有効な組み合わせと無効な組み合わせを個別に計算することでした。残念ながら、合計エラーのため、次のステートメントは true です。

calc_non_skipped() != calc_all() - calc_skipped()

初期値から目的の結果を得る便利で既知の高性能な方法はありますか?

この質問はかなり複雑で、おそらく関連性が限られていることを私は知っています。それにもかかわらず、私の問題を解決するのにいくつかの明快な答えが役立つことを願っています.

現在の実装

現在、これは OpenMP を使用した CPU コードとして実装されています。combo最初に、計算が必要なすべての P を格納する上記の s のベクトルを設定し、それを並列 for ループに渡します。各スレッドにはプライベート結果ベクトルが提供され、並列領域の最後にあるクリティカル セクションが適切な合計に使用されます。

4

2 に答える 2

6

最初に、なぜ N=7 に対して 27 が得られたのか、しばらく戸惑いました(N**2 - N)/2...しかし、インデックス 0 ～ 7、N=8、および P には28の要素があります。日。:-)

しかし、解決策の可能性について:他の目的のために配列 P を保持する必要がありますか? そうでない場合は、それぞれの長さが N の 2 つの中間配列だけで、必要な結果を得ることができると思います。1 つは行の合計用で、もう 1 つは列の合計用です。

これは、あなたがやろうとしていること（ subroutine direct_approach()）と、中間配列を使用して同じ結果を得る方法（ subroutine ）の簡単な例ですrefined_approach()。

#include <cstdlib>
#include <cstdio>

const int N = 7;
const float input_values[N] = { 3.0F, 5.0F, 7.0F, 11.0F, 13.0F, 17.0F, 23.0F };
float P[N][N];      // Yes, I'm wasting half the array.  This way I don't have to fuss with mapping the indices.
float result1[N] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };
float result2[N] = { 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F, 0.0F };

float f(float arg1, float arg2)
{
    // Arbitrary computation
    return (arg1 * arg2);
}

float compute_result(int index)
{
    float row_sum = 0.0F;
    float col_sum = 0.0F;
    int row;
    int col;

    // Compute the row sum
    for (col = (index + 1); col < N; col++)
    {
        row_sum += P[index][col];
    }

    // Compute the column sum
    for (row = 0; row < index; row++)
    {
        col_sum += P[row][index];
    }

    return (row_sum - col_sum);
}

void direct_approach()
{
    int row;
    int col;

    for (row = 0; row < N; row++)
    {
        for (col = (row + 1); col < N; col++)
        {
            P[row][col] = f(input_values[row], input_values[col]);
        }
    }

    int index;

    for (index = 0; index < N; index++)
    {
        result1[index] = compute_result(index);
    }
}

void refined_approach()
{
    float row_sums[N];
    float col_sums[N];
    int index;

    // Initialize intermediate arrays
    for (index = 0; index < N; index++)
    {
        row_sums[index] = 0.0F;
        col_sums[index] = 0.0F;
    }

    // Compute the row and column sums
    // This can be parallelized by computing row and column sums
    //  independently, instead of in nested loops.
    int row;
    int col;

    for (row = 0; row < N; row++)
    {
        for (col = (row + 1); col < N; col++)
        {
            float computed = f(input_values[row], input_values[col]);
            row_sums[row] += computed;
            col_sums[col] += computed;
        }
    }

    // Compute the result
    for (index = 0; index < N; index++)
    {
        result2[index] = row_sums[index] - col_sums[index];
    }
}

void print_result(int n, float * result)
{
    int index;

    for (index = 0; index < n; index++)
    {
        printf("  [%d]=%f\n", index, result[index]);
    }
}

int main(int argc, char * * argv)
{
    printf("Data reduction test\n");

    direct_approach();

    printf("Result 1:\n");
    print_result(N, result1);

    refined_approach();

    printf("Result 2:\n");
    print_result(N, result2);

    return (0);
}

各中間値はほとんどの入力の関数であるため、計算の並列化はそれほど簡単ではありません。合計を個別に計算できますが、それは f(...) を複数回実行することを意味します。N の値が非常に大きい場合に考えられる最良の提案は、より多くの中間配列を使用し、結果のサブセットを計算してから、部分配列を合計して最終的な合計を算出することです。疲れていないときは、そのことを考えなければなりません。

スキップの問題に対処するには: 「入力値 x、y、および z を使用しない」という単純な問題であれば、x、y、および z を do_not_use 配列に格納し、ループして計算するときにそれらの値をチェックできます。合計。スキップする値が行と列の関数である場合、それらをペアとして保存し、ペアを確認できます。

これが解決策のアイデアになることを願っています!

目が覚めたので更新します。 「スキップ」の処理は、スキップする必要があるデータに大きく依存します。最初のケースの別の可能性 - 「入力値 x、y、および z を使用しない」 - 大規模なデータ セットのはるかに高速なソリューションは、間接的なレベルを追加することです。さらに別の配列を作成します。これは整数インデックスの 1 つです。適切な入力のインデックスのみを保存します。たとえば、無効なデータが入力 2 と 5 にある場合、有効な配列は次のようになります。

int valid_indices[] = { 0, 1, 3, 4, 6 };

array を Interrate し、valid_indicesこれらのインデックスを使用して入力配列からデータを取得し、結果を計算します。もう一方の足では、スキップする値が P 配列の両方のインデックスに依存している場合、何らかのルックアップを回避する方法がわかりません。

並列化に戻る - 何があっても、(N**2 - N)/2 回の f() の計算を処理することになります。1 つの可能性は、f() の計算が 2 つの加算よりもかなり長い時間がかかる場合、合計配列の競合が発生することを受け入れることです。非常に多数の並列パスに到達すると、競合が再び問題になりますが、並列パスの数と f() の計算に必要な時間とのバランスを取る「スイート スポット」が存在するはずです。

競合がまだ問題である場合は、いくつかの方法で問題を分割できます。1 つの方法は、行または列を一度に計算することです。行ごとに、各列の合計を個別に計算し、行の合計ごとに現在の合計を保持できます。

もう 1 つのアプローチは、データ空間を分割し、計算をサブセットに分割することです。各サブセットには、独自の行と列の合計配列があります。各ブロックが計算された後、独立した配列を合計して、結果の計算に必要な値を生成できます。

于 2013-06-27T03:49:36.980 に答える
3

これはおそらく素朴で役に立たない答えの1つですが、役立つかもしれません. 私がまったく完全に間違っており、全体のことを誤解していることを遠慮なく言ってください。

では、どうぞ！

基本的な問題

結果関数を少し異なる方法で定義でき、中間値から少なくともいくつかの競合を取り除くことができるように思えます。あなたのP行列が下三角であるとしましょう。上の三角形を下の値の負の値で (仮想的に) 埋めた場合 (および主対角線をすべてゼロで埋めた場合)、結果の各要素を 1 つの行の合計として再定義できます。 、および where-iは、 としてマークされたセル内の値の負の値を意味しますi)

 P     0    1    2    3 
   |--------------------
  0|   x   -0   -1   -3
   |
  1|   0    x   -2   -4
   |
  2|   1    2    x   -5
   |
  3|   3    4    5    x

この行列の各行の合計を計算するために独立したスレッド (同じカーネルを実行) を起動すると、各スレッドは単一の結果要素を書き込みます。問題のサイズが大きすぎて、ハードウェア スレッドが飽和状態になり、ビジーなままになっているようです。

f(x, y)もちろん、注意点は、それぞれを 2 回計算することです。それがどれくらいの費用がかかるのか、またはメモリの競合によって以前にどれだけの費用がかかっていたのかはわかりません。そのため、これが価値のあるトレードオフであるかどうかを判断することはできません。しかしf、本当に高価でない限り、そうかもしれません。

値をスキップする

P計算で無視する必要があるマトリックスの要素が何万もある可能性があると言及しています（事実上それらをスキップします）。

上記で提案したスキームを使用するには、スキップされた要素を(row, col)ペアとして保存する必要があり、各座標ペアの転置も追加する必要があると思います (したがって、スキップされた値の数が 2 倍になります)。のスキップリストの例になりP[6], P[14] and P[18]ます。P(4,0), P(5,4), P(6,3)P(4,0), P(5,4), P(6,3), P(0,4), P(4,5), P(3,6)

次に、このリストを最初に行に基づいて並べ替え、次に列に基づいて並べ替えます。これにより、リストは になりますP(0,4), P(3,6), P(4,0), P(4,5), P(5,4), P(6,3)。

仮想行列の各行がP1 つのスレッド (またはカーネルの単一インスタンスなど) によって処理される場合、スキップする必要がある値を渡すことができます。個人的には、これらすべてを大きな 1D 配列に格納し、各スレッドが参照する必要がある最初と最後のインデックスを渡すだけです (また、渡した最終配列に行インデックスを格納しません。上記の例では、N = 8 の場合、各スレッドに渡される begin と end のペアは次のようになります。 STL のように、空のリストは begin == end で示されます)

Thread 0: 0..1
Thread 1: 1..1 (or 0..0 or whatever)
Thread 2: 1..1
Thread 3: 1..2
Thread 4: 2..4
Thread 5: 4..5
Thread 6: 5..6
Thread 7: 6..6

ここで、各スレッドは続けて、行内のすべての中間値を計算して合計します。列のインデックスをステップ実行している間、スキップされた値のこのリストをステップ実行し、リストに表示される列番号をスキップしています。これは明らかに効率的で簡単な操作です (リストも列ごとにソートされるため、マージのようなものです)。

疑似実装

私は CUDA については知りませんが、OpenCL を使用した経験はあります。インターフェイスは似ていると思います (対象となるハードウェアが同じであるため)。行の処理を行うカーネルの実装を次に示します (つまり、result疑似 C++で ) の 1 つのエントリを計算します。

double calc_one_result (
    unsigned my_id, unsigned N, double const init_values [],
    unsigned skip_indices [], unsigned skip_begin, unsigned skip_end
)
{
    double res = 0;
    for (unsigned col = 0; col < my_id; ++col)
        // "f" seems to take init_values[column] as its first arg
        res += f (init_values[col], init_values[my_id]);
    for (unsigned row = my_id + 1; row < N; ++row)
        res -= f (init_values[my_id], init_values[row]);
    // At this point, "res" is holding "result[my_id]",
    // including the values that should have been skipped
    
    unsigned i = skip_begin;
    // The second condition is to check whether we have reached the
    // middle of the virtual matrix or not
    for (; i < skip_end && skip_indices[i] < my_id; ++i)
    {
        unsigned col = skip_indices[i];
        res -= f (init_values[col], init_values[my_id]);
    }
    for (; i < skip_end; ++i)
    {
        unsigned row = skip_indices[i];
        res += f (init_values[my_id], init_values[row]);
    }
    
    return res;
}

次の点に注意してください。

  1. init_valuesと 関数のセマンティクスはf、質問で説明されているとおりです。

  2. この関数は、result配列内の 1 つのエントリを計算します。具体的には、 を計算するので、これのインスタンスをresult[my_id]起動する必要があります。N

  3. それが書き込む唯一の共有変数はresult[my_id]. 上記の関数は何も書き込みませんが、CUDA に変換すると、最後にそれを書き込む必要があると思います。ただし、 の特定の要素に書き込む人は他にいresultないため、この書き込みによってデータ競合の競合が発生することはありません。

  4. 2 つの入力配列は、init_valuesこのskipped_indices関数の実行中のすべてのインスタンス間で共有されます。

  5. データへのすべてのアクセスは、避けられないと私が信じているスキップされた値を除いて、線形かつシーケンシャルです。

  6. skipped_indices各行でスキップする必要があるインデックスのリストが含まれています。その内容と構造は上記のとおりですが、1 つの小さな最適化が行われています。必要がなかったので、行番号を削除し、列のみを残しました。行番号はとにかく関数に渡され、各呼び出しで使用される配列my_idのスライスは andを使用して決定されます。skipped_indicesskip_beginskip_end

    上記の例では、 のすべての呼び出しに渡される配列は次のcalc_one_resultようになります[4, 6, 0, 5, 4, 3]。

  7. ご覧のとおり、ループを除いて、このコードの唯一の条件分岐はskip_indices[i] < my_id3 番目の for ループにあります。これは無害で完全に予測可能だと思いますが、この分岐でさえコード内で簡単に回避できます。スキップされたアイテムが主対角線を横切る場所を示す別のパラメータを渡す必要があるだけですskip_middle(つまり、行 #my_idの場合、インデックス atskipped_indices[skip_middle]は より大きい最初のものmy_idです)。

結論は

私は決して CUDA と HPC の専門家ではありません。しかし、あなたの問題を正しく理解していれば、この方法でメモリの競合をすべて解消できると思います。また、これが（それ以上の）数値安定性の問題を引き起こすとは思いません。

これを実装するためのコストは次のとおりです。

  • 合計で 2 倍の回数を呼び出しfます (また、いつ呼び出されたかを追跡しrow < colて、結果に を掛けることができ-1ます)。
  • スキップされた値のリストに 2 倍の項目を格納します。このリストのサイズは数千 (数十億ではありません!) であるため、大きな問題にはなりません。
  • スキップされた値のリストを並べ替えます。これもサイズが大きいため問題ありません。

（更新：疑似実装セクションを追加しました。）

于 2013-07-01T07:37:23.030 に答える