OMPスレッド数を増やすと出力値が変更されます

Question

OMPを使用して並列化すると異なる出力が得られるコードを添付します。 これは起こるつもりではありません。

私は、通常の落とし穴を「私的な」変数と「削減」節に適用したと思います。

関連する機能をグローバルコードから添付します。どんな知恵にも感謝しています。

EDIT * MatIplusIplusのi番目の反復とMatPtJplusJplusのj番目の反復から取り出された列ベクトルのarma :: cross productに問題があると思われますか？ *あなたはすべての行列に多くの競合状態を持っている

double calcul(arma::mat MatPt, int iS, int THREADS){ 

int i,j, count = 1; 
double Wr = 0 , Omega; 

arma::mat MatJplusIplus (3 , iS) ; 
arma::mat MatJI (3 , iS) ; 
arma::mat MatJplusI (3 , iS) ; 
arma::mat MatJIplus (3 , iS) ; 

arma::mat MatIplusI (3 , iS) ; 
arma::mat MatJplusJ (3 , iS) ; 

#pragma omp parallel for private(j) reduction(+: Wr) reduction(*: Omega) 
for(i = 0 ; i < iS- (3) ; i++){//rethink boundaries. 

    MatIplusI (0 , i) = MatPt (0, i+1) - MatPt (0, i) ; 
    MatIplusI (1 , i) = MatPt (1, i+1) - MatPt (1, i) ; 
    MatIplusI (2 , i) = MatPt (2, i+1) - MatPt (2, i) ; 

    for(j= i+2 ; j < iS-(1) ; j++){ 

     MatJIplus (0, j) = MatPt (0, j) - MatPt (0, i+1) ; 
     MatJIplus (1, j) = MatPt (1, j) - MatPt (1, i+1) ; 
     MatJIplus (2, j) = MatPt (2, j) - MatPt (2, i+1) ; 

     MatJplusIplus (0, j) = MatPt (0, j+1) - MatPt (0, i+1) ; 
     MatJplusIplus (1, j) = MatPt (1, j+1) - MatPt (1, i+1) ; 
     MatJplusIplus (2, j) = MatPt (2, j+1) - MatPt (2, i+1) ; 

     MatJI (0, j) = MatPt (0, j) - MatPt (0, i) ; 
     MatJI (1, j) = MatPt (1, j) - MatPt (1, i) ; 
     MatJI (2, j) = MatPt (2, j) - MatPt (2, i) ; 

     MatJplusI (0, j) = MatPt (0, j+1) - MatPt (0, i) ; 
     MatJplusI (1, j) = MatPt (1, j+1) - MatPt (1, i) ; 
     MatJplusI (2, j) = MatPt (2, j+1) - MatPt (2, i) ; 

     arma::vec n1 = arma::cross(MatJI.col(j) , MatJplusI.col(j)) ; 
     arma::vec n2 = arma::cross(MatJplusI.col(j) , MatJplusIplus.col(j)) ; 
     arma::vec n3 = arma::cross(MatJplusIplus.col(j) , MatJIplus.col(j)) ; 
     arma::vec n4 = arma::cross(MatJIplus.col(j) , MatJI.col(j)) ; 

     //normalise vectors 
     arma::vec N1 = normalise(n1) ; 
     arma::vec N2 = normalise(n2) ; 
     arma::vec N3 = normalise(n3) ; 
     arma::vec N4 = normalise(n4) ; 

     //take dot product 
     const double ndot1 = arma:: dot (N1 , N2) ; 
     const double ndot2 = arma:: dot (N2 , N3) ; 
     const double ndot3 = arma:: dot (N3 , N4) ; 
     const double ndot4 = arma:: dot (N4 , N1) ; 

     Omega = asin(ndot1) + asin(ndot2) + asin(ndot3) + asin(ndot4) ; 

     MatJplusJ(0 , j) = MatPt (0, j+1) - MatPt (0, j) ; 
     MatJplusJ (1 , j) = MatPt (1, j+1) - MatPt (1, j) ; 
     MatJplusJ (2 , j) = MatPt (2, j+1) - MatPt (2, j) ; 

     arma::vec v = arma::cross(MatJplusJ.col(j) , MatIplusI.col(i)); 

     const double wij = arma:: dot(v, MatJI.col(j)) ; 

     if (wij < 0){ 

      Omega *= -1/ (4* pi) ; 

     } 

     else{ 

      Omega *= 1/ (4* pi) ; 

     } 

     if (Omega != Omega){//incase something goes wrong, ignore ! 

      Omega = 0 ; 

     } 

     Wr += Omega ; 

     //cout << i << " " << j << endl ; 

     //PrivEnd 

    } 

} 

Wr *= 2 ; 


return Wr ; 

}

感謝

Answer 1

はjにアクセスします。例えば、MatJIplus (0, j) = ...は、外側ループの異なる反復によって同時に設定されます。

これらのデータ構造は間接的であるため、それらをプライベートとして宣言するだけでは、ワーカースレッド自体でそれらを構築する必要がありません。各ループ本体内に記述されている行列の部分だけを使用するので、ループ本体の中に宣言することができます。しかし、すべてのループで何度も何度もそれらを再割り当てを避けるために、次のようなparallel forを分割することができます：

arma::mat MatIplusI (3 , iS) ; 
#pragma omp parallel 
{ 
    // Matrices will be implicitly private and correctly constructed locally for each thread. 
    arma::mat MatJplusIplus (3 , iS) ; 
    arma::mat MatJI (3 , iS) ; 
    arma::mat MatJplusI (3 , iS) ; 
    arma::mat MatJIplus (3 , iS) ; 
    arma::mat MatJplusJ (3 , iS) ; 
    #pragma omp parallel for private(j) reduction(+: Wr) private(Omega) 
    { 
     for(i = 0 ; i < iS- (3) ; i++) 

注Omegaの減少は全く理にかなっていないこと、それは代わりに、privateである必要があります。

また、各繰り返し内で列の列のみを使用すると、完全な行列にメモリを割り当てるのはほとんど意味がありません。最初にcolumn(i/j)の内容を保持するのに適したプライベートデータ構造を作成するだけです。