Некоторая ошибка в параллельной реализации Red / Black Successive Over Relaxation (SOR)

Question

Одно из моих назначений - реализовать алгоритм RED / BLACK SOR с использованием MPI.Сетка делится на шахматную доску, и каждая итерация разбивается на две фазы: красную и черную.Во время каждой фазы алгоритм вычисляет либо красные, либо черные неограничные точки сетки.Остальная часть реализации аналогична описанной в wiki .

Полный код: последовательный здесь , параллельный здесь .

Вот последовательная реализация

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

Для параллельной реализацииСетка сначала делится поровну для разных узлов (по столбцам).Например, если размер сетки равен 8x8, а узлов 3, то мы делим сетку как 8x3, 8x3, 8x2 на эти узлы.Во время обмена данными данные обмениваются между левыми и правыми соседями узла.Рисунок ниже должен дать четкую картину всего процесса.

/* 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();

На рисунке показано, как распределена сетка и как соседи общаются.

Проблема в том, что конечный результат параллельного и последовательного SOR, кажется, изменяется на несколько битов.Нужна свежая пара глаз, чтобы просмотреть код для отслеживания этой ошибки, я думаю, что связь между узлами работает нормально.

$ 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

Это выглядит подозрительно:

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

Я бы ожидал что-то вроде:

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

И, возможно, средние два wait () не нужны вообще.

Answer 2

Думали ли вы попробовать на нем отладчик? Я очень предвзят, поскольку я работаю в компании-отладчике MPI, но я только что скачал ваш код и попробовал его в Allinea DDT, и он остановился в stencil (), когда обнаружил чтение за пределами вашего конца. массив, изначально процессом 7.

Чтобы воспроизвести это, скомпилируйте свой код с помощью «-g» для поддержки отладки, и получите Allinea DDT от http://www.allinea.com,, установите его и включите отладку памяти (и убедитесь, что у вас настроены «Защитные страницы») Выше, с 1 страницей в настройках отладки памяти).

Запустите вашу программу, скажем, с 8 процессами, и вы скоро найдете ответ, в считанные секунды.

Удачи!

David