0

コード：

double x(){return (double)rand()/(double)RAND_MAX;}
double y(){return (double)rand()/(double)RAND_MAX;}
double z(){return (double)rand()/(double)RAND_MAX;}

int d(double x, double y, double z){
        if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
        return 0;
    }

double f(double x, double y, double z){
        return 1;
    }




#pragma omp parallel default(none) private(id,numt,j,local_sum,local_good_dots,local_coi,x_,y_,z_) shared(total_sum,good_dots,count_of_iterations)
    {
        local_coi = count_of_iterations;
        id = omp_get_thread_num() + 1;
        numt = omp_get_num_threads();
        #pragma omp for
        for (j = 1; j <= local_coi;  j++){
            x_=x();
            y_=y();
            z_=z();
            if (d(x_,y_,z_) == 1){
                local_sum += f(x_,y_,z_);
                local_good_dots += 1;

            }
        }

        #pragma omp critical
        {
            total_sum = total_sum + local_sum;
            good_dots = good_dots + local_good_dots;
        }
    }

f()コメント：このコードは、面積の関数の3次元積分を計算するためのモンテカルロ法を実現したものですd()。

このコードはマルチスレッドモード（openmp）でより高速に動作することを期待しています。

しかし、何かがうまくいかない。

数時間の変更（reductionopenmpプラグマ、if-conditionの簡略化（のようなf(x_,y_,z_) * d(x_,y_,z_)））の後、私は理解していませんでした。なぜこの単純なループは、スレッドの数が多いほど遅くなるのですか。

しかし、ループの前に座標ごとに3次元配列を生成してドロップするとshared、プログラムはより高速になります。

だから、質問：

このコードを変更する方法と、並列ブロックで許可されている関数（操作）はどれですか？

PS：私が見るように、そのrand機能は許可されていません（または私は間違っていますか？）

手伝ってくれてありがとう！

変更（@HristoIlievの助けを借りて）

double x(){return (double)rand()/(double)RAND_MAX;}
double y(){return (double)rand()/(double)RAND_MAX;}
double z(){return (double)rand()/(double)RAND_MAX;}

int d(double x, double y, double z){
        if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
        return 0;
    }

double f(double x, double y, double z){
        return 1;
    }


#pragma omp parallel default(none) private(j,local_coi,x_,y_,z_) shared(count_of_iterations) reduction(+:total_sum,good_dots)
    {
        local_coi = count_of_iterations;
        #pragma omp for(prng)
        for (j = 1; j <= local_coi;  j++){                    
        #pragma omp critical(prng)
        {
                x_=x();
                y_=y();
                z_=z();
        }   
            if (d(x_,y_,z_) == 1){
                total_sum += f(x_,y_,z_);
                good_dots += 1;

            }
        }
    }
4

1 に答える 1

2

乱数ジェネレーターrand()は、静的に割り当てられたグローバルな状態を使用し、すべてのスレッドで共有されるため、スレッドセーフではありません。複数のスレッドからそれを使用すると、共有変数への保護されていないアクセスという非常に悪いケースに遭遇し、キャッシュを破棄してプログラムの速度を低下させます。使用するrand_r()か、erand48()代わりに使用する必要があります。これらは、提供する必要のある個別の状態ストレージを使用します。private基本的に、スレッドごとに異なるPRNGを作成して、スレッドごとに1つの状態を宣言する必要があります（たとえば、それを持っている）。次に、それに応じてシードする必要があります。そうしないと、統計的に悪い結果が得られます。原則として、1つのジェネレーターの出力を使用してrand48()、他のジェネレーターにシードを設定できます。適度に長い無相関シーケンスを取得するには、十分なはずです。

を使用した実装例を次に示しますrand_r()（これはモンテカルロシミュレーションにとって非常に悪いerand48ジェネレーターであるというわけではありませんが、GNU Scientific Libraryの「MersenneTwister」タイプのジェネレーターを使用するのが最適です）。

unsigned int prng_state;
#pragma omp threadprivate(prng_state)

double x(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double y(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double z(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}

int d(double x, double y, double z){
    if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
    return 0;
}

double f(double x, double y, double z){
    return 1;
}

...

#pragma omp parallel default(none) \
            private(id,numt,x_,y_,z_) \
            shared(count_of_iterations) \
            reduction(+:total_sum,good_dots)
{
    id = omp_get_thread_num() + 1;
    numt = omp_get_num_threads();

    // Sample PRNG seeding code - DO NOT USE IN PRODUCTION CODE!
    prng_state = 67894 + 1337*id;

    #pragma omp for
    for (j = 1; j <= count_of_iterations;  j++){
        x_=x();
        y_=y();
        z_=z();
        if (d(x_,y_,z_) == 1){
            total_sum += f(x_,y_,z_);
            good_dots += 1;
        }
    }
}

これは（品質の観点から）非常に悪い実装ですが、物事がどのように機能するかについてのアイデアを提供するはずです。また、元のコードへの最小限の変更でスレッドセーフを実現する方法でもあります。基本的なポイントは次のとおりです。

  • PRNG状態prng_stateは、OpenMPthreadprivateディレクティブによって各スレッドに対してプライベートになります。
  • rand_r()スレッド固有の状態変数を使用して、、、および;の代わりに使用されますrand()。x()y()z()
  • PRNG状態は、スレッドに依存する方法で初期化されます。たとえばprng_state = 67894 + 1337*id;、異なるスレッドが（うまくいけば）疑似乱数の無相関ストリームを取得するようになります。

rand()とrand_r()は品質が非常に悪いため、これは単なる学術的な例であることに注意してください。PRNGシーケンスが長くなると、さまざまなスレッドで相関ストリームが取得され、統計が損なわれます。を使用してコードを書き直すのはあなたに任せますerand48()。

最初の質問に答えるために-すべてのスレッドセーフな関数呼び出しはparallelブロック内で許可されます。スレッドセーフではない関数を呼び出すこともできますが、（名前付き）構造内の呼び出しを保護する必要がありますcritical。例：

#pragma omp for
for (j = 1; j <= local_coi; j++) {
    #pragma omp critical(prng)
    {
        x_=x();
        y_=y();
        z_=z();
    }
    if (d(x_,y_,z_) == 1) {
        local_sum += f(x_,y_,z_);
        local_good_dots += 1;
    }
}

rand()これにより、への呼び出しが並行して行われないことが保証されます。ただし、共有状態への読み取り-変更-書き込みのようなアクセスが可能であるため、キャッシュ関連の速度が低下します。

また、OpenMPまたは同様の構造を再実装しようとしないでください。reductionコンパイラベンダーは、可能な限り最良の（最速の）方法で実装されるようにするために、すでに多大な努力を払っています。

于 2012-11-19T17:08:14.860 に答える