15

最新の (SSE2+) x86 プロセッサで最大 4096 個の 32 ビット浮動小数点数の配列のソートされたサブセットをマージする高速な方法は何ですか?

以下を想定してください。

  • セット全体のサイズは最大 4096 項目です
  • サブセットのサイズについては議論の余地がありますが、最初は 16 ～ 256 の間であると仮定します。
  • マージで使用されるすべてのデータは、できれば L1 に収まるようにする必要があります
  • L1 データ キャッシュ サイズは 32K です。16K はデータ自体に既に使用されているため、16K で遊ぶことができます。
  • すべてのデータはすでに L1 にあります (可能な限り高い信頼度で) - ソートによって操作されたばかりです
  • すべてのデータは 16 バイトでアラインされています
  • 分岐を最小限に抑えたい (明らかな理由から)

実現可能性の主な基準: In-L1 LSD 基数ソートよりも高速。

上記のパラメータを考慮して、誰かがこれを行うための合理的な方法を知っているかどうかを知りたいです! :)

4

6 に答える 6

8

これは非常に素朴な方法です。(午前 4 時のせん妄に起因する疑似コードのバグはご容赦ください ;)

//4x sorted subsets
data[4][4] = {
  {3, 4, 5, INF},
  {2, 7, 8, INF},
  {1, 4, 4, INF},
  {5, 8, 9, INF}
}

data_offset[4] = {0, 0, 0, 0}

n = 4*3

for(i=0, i<n, i++):
  sub = 0
  sub = 1 * (data[sub][data_offset[sub]] > data[1][data_offset[1]])
  sub = 2 * (data[sub][data_offset[sub]] > data[2][data_offset[2]])
  sub = 3 * (data[sub][data_offset[sub]] > data[3][data_offset[3]])

  out[i] = data[sub][data_offset[sub]]
  data_offset[sub]++


編集:
AVX2 とその収集サポートにより、一度に最大 8 つのサブセットを比較できました。


編集 2:
タイプ キャストによっては、Nehalem での反復ごとに 3 余分なクロック サイクルを削減できる可能性があります (mul: 5、shift+sub: 4)

//Assuming 'sub' is uint32_t
sub = ... << ((data[sub][data_offset[sub]] > data[...][data_offset[...]]) - 1)


編集 3: 2 つ以上の値
を使用することで、特に K が大きくなるにつれて、順不同の実行をある程度悪用できる可能性があります。max

max1 = 0
max2 = 1
max1 = 2 * (data[max1][data_offset[max1]] > data[2][data_offset[2]])
max2 = 3 * (data[max2][data_offset[max2]] > data[3][data_offset[3]])
...
max1 = 6 * (data[max1][data_offset[max1]] > data[6][data_offset[6]])
max2 = 7 * (data[max2][data_offset[max2]] > data[7][data_offset[7]])

q = data[max1][data_offset[max1]] < data[max2][data_offset[max2]]

sub = max1*q + ((~max2)&1)*q


編集4：

コンパイラのインテリジェンスに応じて、三項演算子を使用して乗算を完全に削除できます。

sub = (data[sub][data_offset[sub]] > data[x][data_offset[x]]) ? x : sub


編集5：

コストのかかる浮動小数点比較を避けるためreinterpret_cast<uint32_t*>()に、整数比較になるため、データを単純化することができます。

もう 1 つの可能性は、型付けされていない SSE レジスタを利用し、整数比較命令を明示的に使用することです。

これ< > ==は、バイナリ レベルで float を解釈するときに演算子が同じ結果を生成するために機能します。


編集6：

値の数を SSE レジスタの数に一致させるためにループを十分に展開すると、比較対象のデータをステージングできます。

反復の最後に、選択した最大値/最小値を含むレジスタを再転送し、それをシフトします。

これには索引付けを少しやり直す必要がありますが、ループに を散らかすよりも効率的であることが証明される場合がありますLEA。

于 2012-07-19T02:45:31.303 に答える
3

これはより研究的なトピックですが、d-wayマージソートを使用してブランチの誤予測を最小限に抑えることについて説明しているこの論文を見つけました。

于 2012-07-18T23:46:19.220 に答える
2

SIMDソートアルゴリズムはすでに詳細に研究されています。論文「マルチコアSIMDCPUアーキテクチャでのソートの効率的な実装」では、説明したこと（およびそれ以上）を実行するための効率的なアルゴリズムについて説明しています。

中心的な考え方は、2つの任意の長さのリストのマージをk個の連続する値のブロックのマージに減らすことができるということです（kは4から16の範囲です）。最初のブロックはz[0] = merge(x[0], y[0]).loです。2番目のブロックを取得するために、残りの要素には、からの要素と、からの要素がmerge(x[0], y[0]).hi含まれていることがわかります。ただし、との両方の要素を含めることはできません。これは、複数の要素を含める必要があるためです。したがって、どちらを追加する必要があるかを確認する必要があります。最初の要素が低いものは、必ず最初に表示されますnxxnyynx+ny == kz[1]x[1]y[1]z[1]nx+nyx[1]y[1]z、したがって、これは最初の要素を比較することによって簡単に行われます。そして、マージするデータがなくなるまで、これを繰り返します。

配列が値で終わると仮定した場合の擬似コード+inf：

a := *x++
b := *y++
while not finished:
    lo,hi := merge(a,b)
    *z++ := lo
    a := hi
    if *x[0] <= *y[0]:
        b := *x++
    else:
        b := *y++

（これがマージの通常のスカラー実装にどれほど似ているかに注意してください）

もちろん、実際の実装では条件付きジャンプは必要ありません。たとえば、条件付きでスワップxしyてxorトリックを使用し、無条件に読み取ることができます*x++。

mergeそれ自体は、バイトニックソートで実装できます。ただし、kが低い場合は、命令間の依存関係が多くなり、レイテンシが高くなります。マージする必要のあるアレイの数に応じて、のレイテンシーがマスクされるように十分に高いkを選択できますmerge。または、これが可能な場合は、複数の双方向マージをインターリーブします。詳細については、論文を参照してください。

編集：以下は、k =4の場合の図です。すべての漸近解析では、 kが固定されていると想定しています。

  • 大きな灰色のボックスは、サイズn = m * kの2つの配列をマージしています（写真では、m = 3）。

    ここに画像の説明を入力してください

    1. サイズkのブロックを操作します。
    2. 「ブロック全体のマージ」ボックスは、最初の要素を比較することにより、2つの配列をブロックごとにマージします。これは線形時間操作であり、データをブロックの残りの部分にストリーミングするため、メモリを消費しません。レイテンシーは「merge4」ブロックのレイテンシーによって制限されるため、パフォーマンスは実際には重要ではありません。
    3. 各「merge4」ボックスは2つのブロックをマージし、下位k要素を出力し、上位k要素を次の「merge4」にフィードします。各「merge4」ボックスは制限された数の操作を実行し、「merge4」の数はnで線形です。
    4. したがって、マージの時間コストはnで線形です。また、「merge4」は8つのシリアル非SIMD比較を実行するよりもレイテンシが低いため、非SIMDマージと比較して大幅に高速化されます。

  • 最後に、双方向マージを拡張して多くの配列をマージするために、大きな灰色のボックスを古典的な分割統治法で配置します。各レベルの複雑さは要素数に比例するため、全体の複雑さはO（n log（n / n0））であり、n0は並べ替えられた配列の初期サイズであり、nは最終配列のサイズです。

    ダイアグラム

于 2012-07-22T12:01:23.400 に答える
2

int配列（高価な）ブランチフリーをマージできます。

typedef unsigned uint;
typedef uint* uint_ptr;

void merge(uint*in1_begin, uint*in1_end, uint*in2_begin, uint*in2_end, uint*out){

  int_ptr in [] = {in1_begin, in2_begin};
  int_ptr in_end [] = {in1_end, in2_end};

  // the loop branch is cheap because it is easy predictable
  while(in[0] != in_end[0] && in[1] != in_end[1]){
    int i = (*in[0] - *in[1]) >> 31;
    *out = *in[i];
    ++out;
    ++in[i];
  }

  // copy the remaining stuff ...
}

(*in[0] - *in[1]) >> 31 は *in[0] - *in[1] < 0 と同等であり、これは *in[0] < *in[1] と同等であることに注意してください。代わりにビットシフトトリックを使用して書き留めた理由

int i = *in[0] < *in[1];

すべてのコンパイラが < バージョンのブランチ フリー コードを生成するわけではないということです。

残念ながら、*in[0] < *in[1] ブランチフリーを実際に実装する方法がわからないため、最初はショーストッパーのように見える int の代わりに float を使用しています。ただし、ほとんどの最新のアーキテクチャでは、正の浮動小数点数 (NAN、INF、またはそのような奇妙なものでもありません) のビットパターンを int として解釈し、< を使用してそれらを比較すると、正しい結果が得られます。おそらく、この観察を任意のフロートに拡張します。

于 2012-07-21T21:51:41.310 に答える
2

頭に浮かぶ最も明白な答えは、ヒープを使用した標準的な N-way マージです。それは O(N log k) になります。サブセットの数は 16 ～ 256 であるため、最悪の場合の動作 (それぞれ 16 アイテムの 256 サブセット) は 8N になります。

キャッシュの動作は、完全ではありませんが、合理的であるべきです。ほとんどのアクションが存在するヒープは、おそらくずっとキャッシュに残ります。書き込まれる出力配列の部分もキャッシュにある可能性が高いです。

16K のデータ (ソートされたサブシーケンスを含む配列)、ヒープ (1K、最悪の場合)、およびソートされた出力配列 (再び 16K) があり、32K のキャッシュに収まるようにしたいとします。問題のように聞こえますが、おそらくそうではありません。スワップ アウトされる可能性が最も高いデータは、挿入ポイントが移動した後の出力配列の先頭です。ソートされたサブシーケンスがかなり均一に分散されていると仮定すると、それらをキャッシュに保持するのに十分な頻度でアクセスする必要があります。

于 2012-07-20T03:15:48.103 に答える
1

K 個のリストをマージする単純なマージ カーネルを実行できます。

float *input[K];
float *output;

while (true) {
  float min = *input[0];
  int min_idx = 0;
  for (int i = 1; i < K; i++) {
    float v = *input[i];
    if (v < min) {
      min = v;     // do with cmov
      min_idx = i; // do with cmov
    }
  }
  if (min == SENTINEL) break;
  *output++ = min;
  input[min_idx]++;
}

ヒープがないので、とてもシンプルです。悪い部分は、それが O(NK) であることです。これは、K が大きい場合に悪い可能性があります (O(N log K) であるヒープ実装とは異なります)。したがって、最大の K を選択し (4 または 8 が適切であり、内側のループをアンロールできます)、カスケード マージによってより大きな K を実行します (リストのグループの 8 方向のマージを実行して K=64 を処理し、次に結果の 8 通りのマージ)。

于 2012-07-21T04:34:16.033 に答える