サイズ N (1 ≤M≤N) の母集団から置換なしで K>>1 の単純なランダム長 M のサンプルを生成する。
問題の定義と詳細は、[SRSWOR] を参照してください。
部分 Fisher-Yates シャッフル・アルゴリズム [KnuthV2] の次の実装と oneMKL 乱数ジェネレーター (RNG) を使用して各サンプルを生成します。
A2.1: (初期化ステップ) PERMUT_BUF は自然数1, 2, ..., N を含むものとする
A2.2: i は1 から M までの繰り返し: A2.3: {i,...,N} で一様なランダムな整数 X を生成
A2.4: PERMUT_BUF[i] とPERMUT_BUF[X] を交換
A2.5: (コピーステップ) i は1 から M までの繰り返し : RESULTS_ARRAY[i]=PERMUT_BUF[i]
終了。アルゴリズムを実装するプログラムは 11 969 664 の実験を行います。各実験 (1 から N までの M の一様なランダム自然数のシーケンスを生成) は、実際には N 要素の母集団から部分長 M をランダムにシャッフルします。アルゴリズムのメインループは実際に抽選を行うため、プログラムで各実験は「N から M の抽選」と呼ばれています。
プログラムは M=6 および N=49 を使用して、結果のサンプル (長さ M のシーケンス) を単一配列 RESULTS_ARRAY に格納し、利用可能な並列スレッドをすべて使用します。
ソースコード: サンプル (https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip (英語)) の lottery6of49 フォルダーを参照してください。
#pragma omp parallel
{
thr = omp_get_thread_num(); /* スレッドのインデックス */
VSLStreamStatePtr stream;
/* このスレッドのRNG ストリームを初期化 */
vslNewStream( &stream, VSL_BRNG_MT2203+thr, seed );
...
/* (スレッド番号 thr で) 実験サンプルを生成 */
vslDeleteStream( &stream );
}コードは、OpenMP* #pragma parallel ディレクティブを使用して、すべての CPU のすべての利用可能なプロセッサー・コアを利用します。実験結果の配列 RESULTS_ARRAY は THREADS_NUM の部分に分割されます。ここで、THREADS_NUM は利用可能なスレッド数で、各スレッド (並列領域) は配列の個別の部分を処理します。
oneMKL 基本乱数ジェネレーターは、VSL_BRNG_MT2203 パラメーターを指定することにより、各スレッドで並列の独立したストリームをサポートします。
/* A2.1: 初期化ステップ */
/* PERMUT_BUF は自然数1, 2, ..., N を含むものとする */
for( i=0; i<N; i++ ) PERMUT_BUF[i]=i+1; /* セット{1,...,N} を使用 */
for( sample_num=0; sample_num<EXPERIM_NUM/THREADS_NUM; sample_num++ ){
/* 次の抽選サンプルを生成 (ステップ A2.2、A2.3、および A2.4): */
Fisher_Yates_shuffle(...);
/* A2.5: コピーステップ */
for(i=0; i<M; i++)
RESULTS_ARRAY[thr*ONE_THR_PORTION_SIZE + sample_num*M + i] = PERMUT_BUF[i];
}このコードは各スレッドで部分 Fisher-Yates シャッフル・アルゴリズムを実装します。
多くの実験をシミュレーションする場合、初期化ステップは各実験の最初に 1 回のみ必要です。(実際の抽選のように) PERMUT_BUF 配列の自然数 1...N の順序は重要ではありません。
/* A2.2: i は 0 から M-1 の繰り返し */
Fisher_Yates_shuffle (...)
{
for(i=0; i<M; i++) {
/* A2.3: {i,...,N-1} からランダム自然数X を生成 */
j = Next_Uniform_Int(...);
/* A2.4: PERMUT_BUF[i] とPERMUT_BUF[X] を交換 */
tmp = PERMUT_BUF[i];
PERMUT_BUF[i] = PERMUT_BUF[j];
PERMUT_BUF[j] = tmp;
}
}ループ A2.2 の各反復は、実際に抽選を行います。残りのアイテム PERMUT_BUF[i], ..., PERMUT_BUF[N] を含む箱からランダムなアイテム X を取り出し、アイテム X を結果行 PERMUT_BUF[1],...,PERMUT_BUF[i] の最後に追加します。長さ N の完全な置換ではなく、一部の長さ M のみを生成するため、アルゴリズムは部分的です。
アルゴリズムで説明した擬似コードと異なり、プログラムはゼロベースの配列を使用します。
ステップ A2.3 で、プログラムは Next_Uniform_Int 関数を呼び出して次のランダム整数 X を生成し、{i, ..., N-1} で一様にします (詳細はソースコードを参照)。oneMKL のベクトル化された RNG の能力を引き出しつつ、ベクトル化のオーバーヘッドを最小限に抑えるため、ジェネレーターは L1 キャッシュに収まるサイズ RNGBUFSIZE のベクトル D_UNIFORM01_BUF を生成する必要があります。各スレッドは、独自のバッファー D_UNIFORM01_BUF およびそのバッファーで最後に使用された乱数の後を指すインデックス D_UNIFORM01_IDX を使用します。Next_Uniform_Int 関数の最初の呼び出しでは (あるいはバッファーの乱数がすべて使用された場合は)、長さ RNGBUFSIZE とインデックス D_UNIFORM01_IDX をゼロにセットして vdRngUniform 関数を呼び出し、乱数バッファーを生成し直します。
vdRngUniform( ... RNGBUFSIZE, D_UNIFORM01_BUF ... );oneMKL は同じ分布の乱数値ジェネレーターのみ提供しますが、ステップ A2.3 では異なる範囲のランダム整数が必要なため、[0;1) で一様に分布された倍精度乱数をバッファーに格納した後、「整数スケーリング」ステップで、これらの倍精度値を必要な整数範囲に収まるように変換します。
number 0 distributed on {0,...,N-1} = 0 + {0,...,N-1}
number 1 distributed on {1,...,N-1} = 1 + {0,...,N-2}
...
number M-1 distributed on {M-1,...,N-1} = M-1 + {0,...,N-M}
(前のM ステップを繰り返す)
number M distributed on: see (0)
number M+1 distributed on: see (1)
...
number 2*M-1 distributed on: see (M-1)
(前のM ステップを再度繰り返す)
...
など
整数スケーリング/* 整数スケーリング・ステップ */
for(i=0;i<RNGBUFSIZE/M;i++)
for(k=0;k<M;k++)
I_RNG_BUF[i*M+k] = k + (unsigned int)(D_UNIFORM01_BUF[i*M+k] * (double)(N-k));ここで RNGBUFSIZE は M の倍数。
このコードに関するパフォーマンスの注意点は [SRSWOR] を参照してください。
タスク |
ルーチン |
|---|---|
RNG ストリームを作成して初期化する |
vslNewStream |
区間 [0;1) に一様に分散された倍精度数を生成する |
vdRngUniform |
RNG ストリームを削除する |
vslDeleteStream |
64 ビット境界でアライメントされたメモリーバッファーを割り当てる |
mkl_malloc |
mkl_malloc で割り当てられたメモリーを解放する |
mkl_free |
製品とパフォーマンス情報 |
|---|
パフォーマンスは、利用、構成、およびその他の要因によって異なります。詳細については、www.Intel.com/PerformanceIndex (英語) をご覧ください。 改訂 #20201201 |