Optimiser la boucle externe avec OpenMP et une réduction

Question

Je me bats un peu avec une fonction. Le calcul est erroné si je tente de paralléliser la boucle extérieure avec un

#pragma omp parallel reduction(+:det). 

Quelqu'un peut-il me montrer comment le résoudre et pourquoi il est un échec?

// template<class T> using vector2D = std::vector<std::vector<T>>; 

float Det(vector2DF &a, int n) 
{ 
    vector2DF m(n - 1, vector1DF(n - 1, 0)); 

    if (n == 1) return a[0][0]; 
    if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1]; 

    float det = 0; 
    for (int i = 0; i < n; i++) 
    { 
    int l = 0; 
#pragma omp parallel for private(l) 
    for (int j = 1; j < n; j++) 
    { 
     l = 0; 
     for (int k = 0; k < n; k++) 
     { 
     if (k == i) continue; 
     m[j - 1][l] = a[j][k]; 
     l++; 
     } 
    } 
    det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1); 
    } 

    return det; 
} 

Answer 1

Si vous paralléliser la boucle extérieure, il y a une condition de course sur cette ligne:

m[j - 1][l] = a[j][k]; 

vous voudrez probablement aussi un parallel for reduction au lieu d'un parallel reduction.

Le problème est que m est partagé, même si cela ne serait pas nécessaire étant donné qu'il est complètement écrasé dans la boucle interne. toujours déclarer des variables aussi localement que possible, cela permet d'éviter des problèmes avec mal variables partagées, par exemple:

float Det(vector2DF &a, int n) 
{ 
    if (n == 1) return a[0][0]; 
    if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1]; 

    float det = 0; 
    #pragma omp parallel reduction(+:det) 
    for (int i = 0; i < n; i++) 
    { 
    vector2DF m(n - 1, vector1DF(n - 1, 0)); 
    for (int j = 1; j < n; j++) 
    { 
     int l = 0; 
     for (int k = 0; k < n; k++) 
     { 
     if (k == i) continue; 
     m[j - 1][l] = a[j][k]; 
     l++; 
     } 
    } 
    det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1); 
    } 
    return det; 
} 

Maintenant, cela est exact, mais comme m peut être coûteux d'affecter, la performance pourrait bénéficier de ne pas le faire dans chaque itération. Cela peut être fait en divisant parallel et for directives en tant que tels:

float Det(vector2DF &a, int n) 
{ 
    if (n == 1) return a[0][0]; 
    if (n == 2) return a[0][0] * a[1][1] - a[1][0] * a[0][1]; 

    float det = 0; 
    #pragma omp parallel reduction(+:det) 
    { 
    vector2DF m(n - 1, vector1DF(n - 1, 0)); 
    #pragma omp parallel for 
    for (int i = 0; i < n; i++) 
    { 
     for (int j = 1; j < n; j++) 
     { 
     int l = 0; 
     for (int k = 0; k < n; k++) 
     { 
      if (k == i) continue; 
      m[j - 1][l] = a[j][k]; 
      l++; 
     } 
     } 
     det += std::pow(-1.0, 1.0 + i + 1.0) * a[0][i] * Det(m, n - 1); 
    } 
    } 
    return det; 
} 

Maintenant, vous pouvez aussi simplement déclarer m comme firstprivate, mais cela suppose que le constructeur de copie fait une copie en profondeur complètement indépendant et donc rendre le code plus difficile à raisonner.

Veuillez noter que vous devez toujours inclure la sortie attendue, la sortie réelle et un minimal complete and verifiable example.