👲🏽 👇🏽 👸🏾 GPUでの並列非並列検索またはプライム検索 😌 🐂 👼🏽

議論の暑さの中、すばらしい夏の夜、私はあなたがCUDAにエラトステネスの高速で動くふるいを書くことができることに気付いて愚かでした。 N = 1,000,000,000（9つのゼロ）をターゲットとして。そして伝説が始まりました...

アルゴリズムの詳細については説明しません。たとえば、ここで読むことができ、その時点で持っていたコードをすぐに表示できます。

#include <iostream> #include <math.h> using namespace std; int main() { double number = 1000000000; bool* a = new bool[int(number/2)]; int i,j,result; for (i=0; i<number/2; i++) a[i] = true; for (i=3; i<=floor(sqrt(number)); i+=2) if (a[i/2]) for (j=i*i; j<=number; j+=i*2) a[j/2]=false; result = 0; for (i=0; i<number/2; i++) if (a[i]) result++; cout << result << endl; delete[] a; return 0; }

Core i3 330Mで14〜15秒間実行され、大量のメモリを消費する、わずかに最適化されたシングルスレッドコード。彼から始めましょう。

ブロックふるい

並列化を成功させる鍵は、タスクを独立した部分に分割して、並列アクセス時に発生するメモリの問題を解消することです。そして、ここでは、いわゆるブロックふるいが非常に便利です。これにより、メモリを節約し、大きなふるいをいくつかの小さなふるいに分割することができます。おそらく、これはふるいの興味深い特性によるものです。つまり、NからNまでのルートからのすべての合成数は、2からNまでのルートまでの素数の倍数になります。各ブロックについてのみ、取り消しを開始する要素を見つけます。

最初に頭に浮かぶのは、除数の列挙であり、除数の場合、除数のステップと交差し始めます。ただし、概算によると、このようなアプローチは、可能な最適化をすべて考慮しても、非常に効果がありません。

N = 10000をとる場合、単純なもの（1229個）を見つけるだけで3万回未満のチェックが必要になります。さらに、各要素が少なくとも1回チェックされることを忘れないでください。

ちなみに、このオプションはCUDAに実装され、GeForce 9600 GTグラフィックスカードをほぼ1時間ハングさせました。

試行番号2

これはすべて終了した可能性があります。しかし、いつも予想外のように、決定が下されました。そして、それは次のように見えました。

N = 100で、ブロックの数が10であり、部品のサイズであるとします。予備ふるいには、2、3、5、7の4つの要素があります。予備ふるいのサイズは部品のサイズと等しいため、2番目のブロックから開始できます。それは私たちに知られている数11で始まり、トリプルは私たちが消す素数の例として役立ちます。彼女の参加で明らかに最適化されたため、デュースは考慮していませんが、それについては後で詳しく説明します。前のブロックは10で終了し、最後の3から9の倍数の要素はブロックの終わりから1つの要素に配置されます。したがって、トリプルからこの「テール」を差し引くと、次のブロックで取り消し線が必要な要素の残りの量（最大12個）を見つけることができます。

そして、予備のふるいからのすべての単純なものについて。さらに、「通常の」ふるいで使用される最適化はここで適用できます。つまり、偶数の要素を削除し、それらを2乗数で削除してメモリを節約します。

理論から実践へ

すべての理論的発展をまとめることが残っています。以下がその結果です。

 for(i=0; i<PARTS; i++) { for(j=0; j<PART_SIZE/2; j++) part[j] = true; for(j=1; j<result1; j++) { offset = 0; if(part1[j]*part1[j]<=(i+1)*PART_SIZE) { offset = i*PART_SIZE+part1[j]-i*PART_SIZE%part1[j]; if(offset%2 == 0) offset += part1[j]; } else break; for(k=offset; k<=(i+1)*PART_SIZE; k+=part1[j]*2) if(part[(k-(i*PART_SIZE)-1)/2]) part[(k-(i*PART_SIZE)-1)/2] = false; } if(i==0) result = j - 1; for(j=0; j<PART_SIZE/2; j++) if(part[j]) result++; }