2015-11-11 215 views
0

也許我的問題的解決方案對於openmp的exprience來說顯而易見,但我沒有。我想使用OpenMP加速下面的子例程:使用OpenMP並行化嵌套循環

void Build_ERIS(vector<double> &eris, vector<Atomic_Orbital> &Basis) 
{ 
    int basis_size = Basis.size(); 
    int m = basis_size*(basis_size+1)/2; 
    eris.resize(m*(m+1)/2); 
    bool compute; 
    std::fill(eris.begin(), eris.end(), 0); 

    int i_orbital,j_orbital, k_orbital,l_orbital, i_primitive, j_primitive, k_primitive,l_primitive,ij,kl, ijkl,ijij,klkl; 
    #pragma omp parallel 
    { 
    #pragma omp for ordered 
    for(i_orbital=0; i_orbital<basis_size; i_orbital++){ 
     for(j_orbital=0; j_orbital<i_orbital+1; j_orbital++){ 
    ij = i_orbital*(i_orbital+1)/2 + j_orbital; 
    for(k_orbital=0; k_orbital<basis_size; k_orbital++){ 
     for(l_orbital=0; l_orbital<k_orbital+1; l_orbital++){ 
     kl = k_orbital*(k_orbital+1)/2 + l_orbital; 
     if (ij >= kl) { 


      ijkl = composite_index(i_orbital,j_orbital,k_orbital,l_orbital); 

      ijij = composite_index(i_orbital,j_orbital,i_orbital,j_orbital); 
      klkl = composite_index(k_orbital,l_orbital,k_orbital,l_orbital); 

      for(i_primitive=0; i_primitive<Basis[i_orbital].contraction.size; i_primitive++) 
     for(j_primitive=0; j_primitive<Basis[j_orbital].contraction.size; j_primitive++) 
      for(k_primitive=0; k_primitive<Basis[k_orbital].contraction.size; k_primitive++) 
      for(l_primitive=0; l_primitive<Basis[l_orbital].contraction.size; l_primitive++) 
       eris[ijkl] += 
      normconst(Basis[i_orbital].contraction.exponent[i_primitive],Basis[i_orbital].angular.l, Basis[i_orbital].angular.m, Basis[i_orbital].angular.n)* 
      normconst(Basis[j_orbital].contraction.exponent[j_primitive],Basis[j_orbital].angular.l, Basis[j_orbital].angular.m, Basis[j_orbital].angular.n)* 
      normconst(Basis[k_orbital].contraction.exponent[k_primitive],Basis[k_orbital].angular.l, Basis[k_orbital].angular.m, Basis[k_orbital].angular.n)* 
      normconst(Basis[l_orbital].contraction.exponent[l_primitive],Basis[l_orbital].angular.l, Basis[l_orbital].angular.m, Basis[l_orbital].angular.n)* 
      Basis[i_orbital].contraction.coef[i_primitive]* 
      Basis[j_orbital].contraction.coef[j_primitive]* 
      Basis[k_orbital].contraction.coef[k_primitive]* 
      Basis[l_orbital].contraction.coef[l_primitive]* 
      ERI_int(Basis[i_orbital].contraction.center.x, Basis[i_orbital].contraction.center.y, Basis[i_orbital].contraction.center.z, Basis[i_orbital].contraction.exponent[i_primitive],Basis[i_orbital].angular.l, Basis[i_orbital].angular.m, Basis[i_orbital].angular.n, 
       Basis[j_orbital].contraction.center.x, Basis[j_orbital].contraction.center.y, Basis[j_orbital].contraction.center.z, Basis[j_orbital].contraction.exponent[j_primitive],Basis[j_orbital].angular.l, Basis[j_orbital].angular.m, Basis[j_orbital].angular.n, 
       Basis[k_orbital].contraction.center.x, Basis[k_orbital].contraction.center.y, Basis[k_orbital].contraction.center.z, Basis[k_orbital].contraction.exponent[k_primitive],Basis[k_orbital].angular.l, Basis[k_orbital].angular.m, Basis[k_orbital].angular.n, 
       Basis[l_orbital].contraction.center.x, Basis[l_orbital].contraction.center.y, Basis[l_orbital].contraction.center.z, Basis[l_orbital].contraction.exponent[l_primitive],Basis[l_orbital].angular.l, Basis[l_orbital].angular.m, Basis[l_orbital].angular.n); 

      /**/ 
     } 
     } 

    } 

     } 
    } 
    } 
} 

我關心的是關於肯定的是,OpenMP並行化後,在鬩神星的減排量的計算[IJKL],仍用相同的值的最佳方式例程的序列版本?我怎樣才能以數字安全的方式進行循環融合?

+0

我不能說我多久看到一個問題題目出現在這裏。我不是[tag:openmp]的專家,但我敢肯定你應該做一點[更多研究](http://stackoverflow.com/search?q=%5Bc%2B%2B%5D% 5Bopenmp%5Dparallelize +循環)確保不要求重複的問題。 –

回答

0

我看到的幾件事。

1)#pragma for ordered表示:按順序執行該循環的每一次迭代。這基本上意味着,當你在「並行」執行時,你的所有工作都將以串行方式完成。去掉它。

2)您還沒有聲明任何變量是共享的還是私有的。請注意,默認情況下,所有變量都將被共享,因此您的案例ijkl可以由任何迭代工作的任何線程訪問。毫無疑問,如果迭代100改變了變量ij,而迭代1認爲它正在使用它,那麼這將會導致競爭狀況。

3)您正確注意的變量eris[ijkl]必須正確減少。如果ijkli_orbital循環中永遠不可能是兩次不同迭代的相同值,那麼您就像原樣;沒有兩個線程可能會同時更改相同的變量eris[ijkl]。如果它可以是相同的值,那麼你have to carefully handle reduction on the array

4)下面是你應該爲初學者工作。這假設ijkl將永遠不會是兩個不同迭代的相同值,並且您的函數不會接受任何非常量引用(可能會將輸入變量假設爲輸出變量)。

#pragma omp parallel for private(i_orbital, j_orbital, ij, k_orbital, l_orbital, kl, ijkl, ijij, klkl, i_primitive, j_primitive, k_primitive, l_primitive)