Un bug dans la mise en œuvre de Redressement successif rouge/noir parallèle (SOR)

Question

Une de mes tâches consiste à implémenter l'algorithme RED/BLACK SOR à l'aide de MPI. La grille est divisée en damier et chaque itération est divisée en deux phases rouge et noire. Au cours de chaque phase, l'algorithme calcule les points non-limites rouges ou noirs de la grille. Le reste de l'implémentation est similaire à comme défini dans wiki. Full Code: séquentiel here, parallèle here.

Voici la mise en oeuvre séquentielle

iteration = 0; 
do { 
    maxdiff = 0.0; 
    for (phase = 0; phase < 2 ; phase++){ 
     for (i = 1 ; i < N-1 ; i++){ 
      for (j = 1 + (even(i)^phase); j < N-1 ; j += 2){ 
       Gnew = stencil(G,i,j); 
       diff = fabs(Gnew - G[i][j]); 
       if (diff > maxdiff) 
        maxdiff = diff; 
       G[i][j] = G[i][j] + omega * (Gnew-G[i][j]); 
      } 
     } 
    } 
    iteration++; 
} while (maxdiff > stopdiff); 

de la grille de mise en œuvre en parallèle est d'abord divisé de manière égale pour les différents noeuds (colonne sage). Par exemple, si la taille de la grille est 8x8 et les nœuds 3, nous divisons la grille en 8x3, 8x3, 8x2 entre ces nœuds. Pendant l'échange, les données sont échangées vers et depuis les voisins gauche et droit du nœud. La figure ci-dessous devrait donner une image claire de l'ensemble du processus.

/* Initializing MPI */ 
MPI_Init(&argc, &argv); 
MPI_Comm_size(MPI_COMM_WORLD, &totalnodes); 
MPI_Comm_rank(MPI_COMM_WORLD, &mynode); 

// divide grid equally among different nodes 
range(1, N - 1, totalnodes, mynode, &jsta, &jend); 

// neigboring nodes 
inext = mynode + 1; 
iprev = mynode - 1; 

iteration = 0; 
do { 
    maxdiff = 0.0; 
    for (phase = 0; phase < 2 ; phase++){ 

     // exchange column with neigboring node 
     exchange(phase); 

     for (i = 1 ; i < N-1 ; i++){   // row 

      for (j = jsta + (even(i)^phase); j <= jend ; j += 2){ // column 

       Gnew = stencil(G,i,j); 
       diff = fabs(Gnew - G[i][j]); 
       if (diff > maxdiff) 
        maxdiff = diff; 
       G[i][j] = G[i][j] + omega * (Gnew-G[i][j]); 
      } 
     } 
    } 

    MPI_Allreduce(&maxdiff, &diff, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD); 
    maxdiff = diff; 

    iteration++; 
} while (maxdiff > stopdiff); 

MPI_Barrier(MPI_COMM_WORLD); 
MPI_Finalize(); 

La figure décrit comment la grille est divisée et les voisins communiquent.

Le problème est que le résultat final du SOR parallèle et séquentiel semble varier de quelques bits. Besoin d'un nouvel oeil de paire pour parcourir le code pour suivre ce bug, je pense que la communication entre les nœuds fonctionnent bien.

$ cc -o sor-seq sor-seq.c -lm  
$ ./sor-seq 8 -print 
    Running 8 x 8 SOR 
    6.006  5.525  5.330  5.251  5.234  5.276  5.417  5.799 
    6.621  6.204  5.984  5.879  5.850  5.892  6.032  6.338 
    6.952  6.687  6.523  6.432  6.395  6.409  6.483  6.640 
    7.181  7.069  6.988  6.931  6.891  6.864  6.852  6.857 
    7.382  7.420  7.429  7.414  7.373  7.306  7.203  7.059 
    7.607  7.799  7.896  7.920  7.884  7.782  7.595  7.294 
    7.926  8.273  8.436  8.488  8.459  8.344  8.101  7.643 
    8.506  8.929  9.088  9.136  9.120  9.033  8.821  8.298 


$ mpicc -o sor-par sor-par.c 
$ mpirun -n 3 ./sor-seq 8 -print 
    Running 8 x 8 SOR 
    5.940  5.383  5.092  4.882  4.677  4.425  4.072  3.507 
    6.496  5.939  5.542  5.201  4.839  4.392  3.794  2.950 
    6.786  6.334  5.938  5.542  5.086  4.512  3.761  2.773 
    6.994  6.672  6.334  5.942  5.450  4.809  3.964  2.873 
    7.197  7.028  6.784  6.442  5.965  5.308  4.414  3.228 
    7.445  7.457  7.335  7.075  6.660  6.045  5.157  3.896 
    7.807  8.020  8.022  7.864  7.555  7.055  6.273  5.032 
    8.443  8.795  8.868  8.805  8.640  8.348  7.848  6.920 

    Node: 0 
     5.940  5.383  5.092  4.882  0.000  0.000  0.000  0.000 
     6.496  5.939  5.542  5.201  0.000  0.000  0.000  0.000 
     6.786  6.334  5.938  5.542  0.000  0.000  0.000  0.000 
     6.994  6.672  6.334  5.942  0.000  0.000  0.000  0.000 
     7.197  7.028  6.784  6.442  0.000  0.000  0.000  0.000 
     7.445  7.457  7.335  7.075  0.000  0.000  0.000  0.000 
     7.807  8.020  8.022  7.864  0.000  0.000  0.000  0.000 
     8.443  8.795  8.868  8.805  0.000  0.000  0.000  0.000 

    Node: 1 
     0.000  0.000  5.092  4.882  4.677  4.425  4.072  0.000 
     0.000  0.000  5.542  5.201  4.839  4.392  3.794  0.000 
     0.000  0.000  5.938  5.542  5.086  4.512  3.761  0.000 
     0.000  0.000  6.334  5.942  5.450  4.809  3.964  0.000 
     0.000  0.000  6.784  6.442  5.965  5.308  4.414  0.000 
     0.000  0.000  7.335  7.075  6.660  6.045  5.157  0.000 
     0.000  0.000  8.022  7.864  7.555  7.055  6.273  0.000 
     0.000  0.000  8.868  8.806  8.640  8.348  7.848  0.000 

    Node: 2 
     0.000  0.000  0.000  0.000  0.000  4.425  4.072  3.507 
     0.000  0.000  0.000  0.000  0.000  4.392  3.794  2.950 
     0.000  0.000  0.000  0.000  0.000  4.512  3.761  2.773 
     0.000  0.000  0.000  0.000  0.000  4.809  3.964  2.873 
     0.000  0.000  0.000  0.000  0.000  5.308  4.414  3.228 
     0.000  0.000  0.000  0.000  0.000  6.045  5.157  3.896 
     0.000  0.000  0.000  0.000  0.000  7.055  6.273  5.032 
     0.000  0.000  0.000  0.000  0.000  8.348  7.848  6.920 

Answer 1

Avez-vous pensé à essayer un débogueur? Je suis très biaisé comme je travaille pour une société de débogage MPI, mais je viens de télécharger votre code, et l'ai essayé dans Allinea DDT, et il a arrêté dans stencil() quand il a détecté une lecture au-delà de la fin de votre tableau Pour reproduire cela, compilez votre code avec "-g" pour la prise en charge du débogage et obtenez Allinea DDT à partir de http://www.allinea.com, installez-le et activez le débogage de la mémoire (et assurez-vous d'avoir configuré "Guard Pages" Ci-dessus, avec 1 page dans les paramètres de débogage de la mémoire). Lancez votre programme, avec, disons, 8 processus, et vous devriez trouver rapidement votre réponse, en quelques secondes.

Bonne chance!

David

Answer 2

Cela semble suspect:

MPI_Isend(bufs1, icnt1, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqs1); 
MPI_Isend(bufs2, icnt2, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqs2); 

MPI_Irecv(bufr1, N, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqr1); 
MPI_Irecv(bufr2, N, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqr2); 

MPI_Wait(&ireqs1, &istatus); 
MPI_Wait(&ireqs2, &istatus); 
MPI_Wait(&ireqr1, &istatus); 
MPI_Wait(&ireqr2, &istatus); 

Je me serais attendu à quelque chose comme:

MPI_Isend(bufs1, icnt1, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqs1); 
MPI_Isend(bufs2, icnt2, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqs2); 

MPI_Wait(&ireqs1, &istatus); 
MPI_Wait(&ireqs2, &istatus); 

MPI_Irecv(bufr1, N, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqr1); 
MPI_Irecv(bufr2, N, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqr2); 

MPI_Wait(&ireqs1, &istatus); 
MPI_Wait(&ireqs2, &istatus); 

Et peut-être au milieu de deux attente (s) ne sont pas nécessaires du tout.