雅各比在MPI放松

我在1小时前创建了一个问题,但问得不好,所以我重新创建了一个问题。

我在C中得到了Jacobi放松的代码:

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

我运行这个:

  • 数组4096×4096
  • iter_max = 1000
  • 错误= 1.0e-6
  • 16核

我已将此代码与OpenACC并行化。 现在,我想使用MPI来了解它是如何工作的。 但是,对于我做的第一个实现,我没有很好的结果(新数组构造不好)。 如何将此代码段与MPI并行化?

这是我为类似案例编写的代码,您可以将其用作指南。

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