Параллельное суммирование с помощью openMP - что делать, если я не могу использовать предложение сокращения?

Question

У меня есть функция внутри класса, которую я называю миллионы раз в моем коде. В этой функции требуются петли, которые могут быть распараллелены. Моя проблема заключается в том, что они выполняют суммы, которые хранятся в переменных, отличных от sclar. вот код.

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; 
      } 
     }    
    }   
} 


} 

Вот пример того, что я сделал в другой части кода:

#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; 
     } 

Спасибо заранее.

Код после @phadjido предложений:

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

Если пользовательские сокращения не поддерживаются компилятором, вы можете просто реализовать операцию сокращения по своему усмотрению. В приведенном ниже коде показано, как это можно сделать для вашего второго примера. Обратите внимание, что pi_2 инициализируется в начале цикла. В вашем примере pi_2 является частной переменной и, возможно, не был инициализирован до нуля. Вам понадобится firstprivate и правильная инициализация pi_2 перед параллельной областью.

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; 
}