Релаксация Якоби в MPI

Я создал вопрос 1 час назад, но он не был правильно задан, поэтому я воссоздал его.

У меня есть код, который является релаксацией Якоби в C:

while ( error > tol && iter < iter_max ) {
    error = 0.0;

    for( int j = 1; j < n-1; j++)
    {   
        for( int i = 1; i < m-1; i++ )
        {   
            Anew[j][i] = 0.25 * ( A[j][i+1] + A[j][i-1]
                                + A[j-1][i] + A[j+1][i]);
            error = fmax( error, fabs(Anew[j][i] - A[j][i]));
        }
    }

    for( int j = 1; j < n-1; j++)
    {   
        for( int i = 1; i < m-1; i++ )
        {   
            A[j][i] = Anew[j][i];
        }
    }

    if(iter % 100 == 0) printf("%5d, %0.6f\n", iter, error);

    iter++;
}

Я запускаю это с помощью:

  • массивы 4096x4096
  • iter_max = 1000
  • ошибка = 1.0e-6
  • 16 ядер

Я распараллелил этот код с OpenACC. Теперь я хочу использовать MPI, чтобы попытаться понять, как это работает. Однако для первых реализаций, которые я сделал, у меня не очень хорошие результаты (новый массив построен плохо). Как я могу распараллелить этот раздел кода с помощью MPI?


c parallel-processing mpi
person Vincent ROSSIGNOL    schedule 22.06.2016    source источник
comment
Если я ввожу термин Релаксация Якоби в MPI в поле поиска широко используемой поисковой системы, первое, что я получаю, — это учебник (с решением) с очень авторитетного сайта, охватывающий простую итерацию Якоби. . Это, вероятно, гораздо лучшая отправная точка для вашей работы, чем задавать этот довольно общий вопрос здесь.   -  person High Performance Mark    schedule 22.06.2016
comment
Я видел этот сайт. Однако я не думаю, что это хорошо работает, если вы хотите создать массив. Он не объединяет все части созданного массива. Кроме того, он предназначен только для 12x12. Но я согласен, что это хорошая отправная точка.   -  person Vincent ROSSIGNOL    schedule 22.06.2016
comment
Добро пожаловать в Stackoverflow! Вас может заинтересовать библиотека PETSc. В частности, функция DMDACreate2d() создает распределенный массив. Например, см. этот пример о неоднородном лапласиане. Распределенные массивы очень распространены в программах MPI. Действительно, он ограничивает объем обмена данными между узлами и обеспечивает хороший баланс между занимаемой памятью и процессорным временем между процессами.   -  person francis    schedule 22.06.2016

Ответы (1)


arrow_upward
0
arrow_downward

Вот код, который я написал для аналогичного случая, и вы можете использовать его в качестве руководства.

do {
  iter++;

  MPI_Irecv(&old[1][0], 1, myHelloVector, nbrs[LEFT], 2, MPI_COMM_WORLD, \
            &requestFourR[LEFT]); 
  MPI_Irecv(&old[1][chunkSize[1]+1], 1, myHelloVector, nbrs[RIGHT], 1, \
            MPI_COMM_WORLD, &requestFourR[RIGHT]); 

  MPI_Irecv(&old[0][1], chunkSize[1], MPI_FLOAT, nbrs[UP], 4, \
            MPI_COMM_WORLD, &requestFourR[UP]);
  MPI_Irecv(&old[chunkSize[0]+1][1], chunkSize[1], MPI_FLOAT, \
            nbrs[DOWN], 3, MPI_COMM_WORLD, &requestFourR[DOWN]);

  MPI_Issend(&old[1][1], 1, myHelloVector, nbrs[LEFT], 1, \
            MPI_COMM_WORLD, &requestFourS[LEFT]);     
  MPI_Issend(&old[1][chunkSize[1]], 1, myHelloVector, nbrs[RIGHT], 2, \
            MPI_COMM_WORLD, &requestFourS[RIGHT]);

  MPI_Issend(&old[1][1], chunkSize[1], MPI_FLOAT, nbrs[UP], 3, \
            MPI_COMM_WORLD, &requestFourS[UP]);
  MPI_Issend(&old[chunkSize[0]][1], chunkSize[1], MPI_FLOAT, nbrs[DOWN], 4, \
            MPI_COMM_WORLD, &requestFourS[DOWN]);

  calImage(old, new, edge, chunkSize[ROWS], chunkSize[COLS]);

  for (itr = 0; itr < 4; itr++) {
    MPI_Waitany(4, &requestFourR[0], &index, &status);
    switch ( index ) {  /* status.MPI_TAG) */
      case 0: /* RIGHT */
              j = 1;
              for (i = 2; i < chunkSize[0]; i++) {
                new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+ \
                          old[i][j+1] - edge[i][j]);
              }
              break;
      case 1: /* LEFT */
              j = chunkSize[1];
              for (i = 2; i < chunkSize[0]; i++) {
                new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+ \
                          old[i][j+1] - edge[i][j]);
              }
              break;
      case 2: /* DOWN */
              i = 1;
              for (j = 2; j < chunkSize[1]; j++) {
                new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+ \
                          old[i][j+1] - edge[i][j]);
              }
              break;
      case 3: /* UP */
              i = chunkSize[0];
              for (j = 2; j < chunkSize[1]; j++) {
                new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+ \
                          old[i][j+1] - edge[i][j]);
              }
              break;
    }
  }

  i = 1; j = 1;
  new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+old[i][j+1] - \
              edge[i][j]);

  i = 1; j = chunkSize[1];
  new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+old[i][j+1] - \
              edge[i][j]);

  i = chunkSize[0]; j = 1;
  new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+old[i][j+1] - \
              edge[i][j]);

  i = chunkSize[0]; j = chunkSize[1];
  new[i][j] = 0.25f*(old[i-1][j]+old[i+1][j]+old[i][j-1]+old[i][j+1] - \
              edge[i][j]);

  MPI_Waitall(4, requestFourS, statusS);

  temp = old;
  old = new;
  new = temp;
} while(your_stopping_condition);

Функция calImage() выполняет вычисления, которые не зависят от операции замены ореола.

void calImage(float **image, float **newImage, float **edge, \
          int rows, int cols) {
  int i, j;

  for (i = 2; i < rows; i++) {
    for (j = 2; j < cols; j++) {
      newImage[i][j] = 0.25f * (image[i-1][j] \
                        + image[i+1][j] \
                        + image[i][j-1] \
                        + image[i][j+1] \
                        - edge[i][j]);
    }
  }
}
person Angelos    schedule 01.07.2016