Parallele Summierung mit openMP - was kann ich tun, wenn ich die Reduktionsklausel nicht nutzen kann?

Question

Ich habe eine Funktion innerhalb einer Klasse, die ich millionenfach in meinem Code aufrufen. In dieser Funktion gibt es anspruchsvolle Schleifen, die parallelisiert werden können. Mein Problem ist, dass sie Summierungen durchführen, die in nicht-skalaren Variablen gespeichert sind. Hier ist der Code.

void Forces::ForRes(vector<vector<double> > & Posicoes,double epsilon,double sigma,double rc,double L) 
{  

double rij_2,rij_6,rij_12,rijx,rijy,rijz; 
double t; 
double sigma6 = pow(sigma,6); 
double sigma12 = pow (sigma6,2);    
for (unsigned int i = 0 ; i < Posicoes.size() - 1 ; i++) 
{   
    for (unsigned int j = i + 1 ; j < Posicoes.size() ; j++) 
    { 
      rijx = (Posicoes[i][0]-Posicoes[j][0]) - L*round((Posicoes[i][0]-Posicoes[j][0])/L); 
      rijy = (Posicoes[i][1]-Posicoes[j][1]) - L*round((Posicoes[i][1]-Posicoes[j][1])/L); 
      rijz = (Posicoes[i][2]-Posicoes[j][2]) - L*round((Posicoes[i][2]-Posicoes[j][2])/L); 
      rij_2 = rijx*rijx + rijy*rijy + rijz*rijz; 
      rij_6 = pow(rij_2,3); 
      rij_12 = pow(rij_6,2); 

     if (rij_2 <= rc*rc) 
     { 
       U += 4*epsilon*((sigma12)/(rij_12)- (sigma6)/(rij_6));         
      for (int k =0 ; k <3 ; k++) 
      {   
       t = ((24*epsilon)/(rij_2))*(2*(sigma12)/(rij_12)- (sigma6)/(rij_6))*((Posicoes[i][k]-Posicoes[j][k]) 
        - L*round((Posicoes[i][k]-Posicoes[j][k])/L)); 

       F[i][k] += t;   
        F[j][k] -= t; 
      } 
     }    
    }   
} 


} 

Hier ist ein Beispiel, das ich in einem anderen Teil des Codes habe:

#pragma omp parallel for default(shared) reduction(+:K) private(pi_2) 
    for (int i = 0 ; i < Nparticulas;i++) 
     { 

     for (int k = 0 ; k < 3 ; k++) 
      { 
       pi_2 += Momentos.M[i][k]*Momentos.M[i][k]; 
      } 

      K += pi_2/2; 
      pi_2 = 0; 
     } 

Vielen Dank im Voraus.

-Code nach @phadjido Vorschlag:

void Forces::ForRes(vector<vector<double> > & Posicoes,double epsilon,double sigma,double rc,double L) 
{ 

    double rij_2,rij_6,rij_12,rijx,rijy,rijz; 
    double t; 
    double sigma6 = pow(sigma,6); 
    double sigma12 = pow (sigma6,2); 

    U = 0; 
    unsigned int j; 

    for (unsigned int i = 0 ; i < Posicoes.size() - 1 ; i++) 
    { 
     #pragma omp parallel private (rij_2,rij_6,rij_12,j) 
     { 

      double Up = 0; 
      vector <vector <double> > Fp(Posicoes.size() , vector<double>(Posicoes[0].size(),0));  

      #pragma omp for 
      for (j = i + 1 ; j < Posicoes.size() ; j++) 
       { 
       rijx = (Posicoes[i][0]-Posicoes[j][0]) - L*round((Posicoes[i][0]-Posicoes[j][0])/L); 
       rijy = (Posicoes[i][1]-Posicoes[j][1]) - L*round((Posicoes[i][1]-Posicoes[j][1])/L); 
       rijz = (Posicoes[i][2]-Posicoes[j][2]) - L*round((Posicoes[i][2]-Posicoes[j][2])/L); 
       rij_2 = rijx*rijx + rijy*rijy + rijz*rijz; 
       rij_6 = pow(rij_2,3); 
       rij_12 = pow(rij_6,2); 

     if (rij_2 <= rc*rc) 
     { 

      Up += 4*epsilon*((sigma12)/(rij_12)- (sigma6)/(rij_6)); 


     for (int k =0 ; k <3 ; k++) 
      {   
       t = ((24*epsilon)/(rij_2))*(2*(sigma12)/(rij_12)- (sigma6)/(rij_6))*((Posicoes[i][k]-Posicoes[j][k]) 
       - L*round((Posicoes[i][k]-Posicoes[j][k])/L)); 
       Fp[i][k] += t; 

       Fp[j][k] -= t; 
       } 
      } 

     } 

     #pragma omp atomic 
     U += Up; 
     for(j = i + 1 ; j < Posicoes.size() ; j++) 
     {  
      for (int k = 0 ; k < 3; k++) 
       {    
      #pragma omp atomic 
      F[i][k] += Fp[i][j]; 
      #pragma omp atomic 
      F[j][k] -= Fp[j][k]; 
       } 

     }   
    } 
    } 

}

Answer 1

Wenn Benutzer definiert Reduktionen werden vom Compiler nicht unterstützt, können Sie einfach die Reduktionsoperation auf eigene Faust umzusetzen. Der folgende Code zeigt, wie dies für Ihr zweites Beispiel getan werden kann. Bitte beachten Sie, dass pi_2 am Anfang der Schleife initialisiert wird. In Ihrem Beispiel ist pi_2 eine private Variable und wurde möglicherweise nicht auf Null initialisiert. Sie müßten zuerst pi_2 vor der parallelen Region initialisieren und richtig initialisieren.

K = 0; 
#pragma omp parallel private(pi_2) 
{ 
    double local_K = 0; /* initialize here */ 

    #pragma omp for 
    for (int i = 0 ; i < Nparticulas;i++) 
    { 
     pi_2 = 0; /* be careful */ 
     for (int k = 0 ; k < 3 ; k++) 
     { 
      pi_2 += Momentos.M[i][k]*Momentos.M[i][k]; 
     } 
     local_K += pi_2/2; 
    } 

    #pragma omp atomic 
     K += local_K; 
}