В вашей программе есть утечка памяти; это:
MPI_Isend(A, Rows, MPI_DOUBLE, my_rank+1, 0, MPI_COMM_WORLD, &request[1]);
MPI_Irecv(AA, Rows, MPI_DOUBLE, my_rank+1, MPI_ANY_TAG, MPI_COMM_WORLD, &request[3]);
MPI_Wait(&request[3], &status[3])
утечка ресурсов, связанных с запросом MPI_Isend. Вы вызываете это Rows*Columns раз за итерацию, по-видимому, за много итераций; но вы звоните только подождите один из запросов. Вероятно, вам нужно выполнить MPI_Waitall() для двух запросов.
Но кроме того, ваша программа очень запутанная. Ни одна разумная MPI-программа не должна иметь такой серии if (rank == ...) операторов. И поскольку вы не выполняете никакой реальной работы между неблокирующими отправкой / получением и ожиданием, я не понимаю, почему вы просто не используете MPI_Sendrecv или что-то в этом роде. Что ваша программа пытается достичь?
UPDATE
Хорошо, похоже, вы делаете стандартную вещь, заполняющую гало. Несколько вещей:
Каждому заданию не нужны свои собственные массивы - A / AA для ранга 0, B / BB для ранга 1 и т. Д. Память распределена, а не разделена; никакие ранги не могут видеть другие массивы, поэтому нет необходимости беспокоиться о перезаписи их. (Если бы было, вам не нужно отправлять сообщения). Кроме того, подумайте, насколько сложнее это работает на разных количествах процессов - вам придется добавлять новые массивы в код каждый раз, когда вы меняете количество используемых процессоров.
Вы можете читать / записывать непосредственно в массив V вместо использования копий, хотя изначально копии могут быть проще для понимания.
Я написал здесь небольшую версию кода, заполняющего гало, используя ваши имена переменных (Tmyo, Nmyo, V, знаки i и y и т. Д.). У каждой задачи есть только часть более широкого массива V, и она обменивается данными своего края только со своими соседями. Он использует символы, чтобы вы могли видеть, что происходит. Он заполняет свою часть массива V своим номером #, а затем обменивается данными своего края со своими соседями.
Я бы СИЛЬНО рекомендовал бы вам сесть с книгой MPI и разобраться с ее примерами. Мне нравится Использование MPI , но есть много других. Есть также много хороших учебников MPI там. Я думаю, что не будет преувеличением сказать, что 95% книг и учебных пособий по MPI (например, наши здесь - см. Части 5 и 6) пройдут именно эту процедуру в качестве одного из своих первых больших проработанных примеров. Они назовут это заполнением гало или заполнением защитной ячейки или обменом границ или чем-то еще, но все это сводится к передаче крайних данных.
#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
char **alloc_2d_char(const int rows, const int cols) {
char *data = (char *)malloc(rows*cols*sizeof(char));
char **array= (char **)malloc(rows*sizeof(char*));
for (int i=0; i<rows; i++)
array[i] = &(data[cols*i]);
return array;
}
void edgeDataFill(char **locV, const int locNmyo, const int locTmyo,
const int ncols, const int myrow, const int mycol,
const int size, const int rank) {
MPI_Datatype leftright, updown;
int left, right, up, down;
int lefttag = 1, righttag = 2;
int uptag = 3, downtag = 4;
MPI_Status status;
/* figure out our neighbours */
left = rank-1;
if (mycol == 0) left = MPI_PROC_NULL;
right = rank+1;
if (mycol == ncols-1) right = MPI_PROC_NULL;
up = rank - ncols;
if (myrow == 0) up = MPI_PROC_NULL;
down = rank + ncols;
if (down >= size) down = MPI_PROC_NULL;
/* create data type for sending/receiving data left/right */
MPI_Type_vector(locNmyo, 1, locTmyo+2, MPI_CHAR, &leftright);
MPI_Type_commit(&leftright);
/* create data type for sending/receiving data up/down */
MPI_Type_contiguous(locTmyo, MPI_CHAR, &updown);
MPI_Type_commit(&updown);
/* Send edge data to our right neighbour, receive from left.
We are sending the edge (locV[1][locTmyo]..locV[locNmyo][locTmyo]),
and receiving into edge (locV[0][1]..locV[locNmyo][locTmyo]) */
MPI_Sendrecv(&(locV[1][locTmyo]), 1, leftright, right, righttag,
&(locV[1][0]), 1, leftright, left, righttag,
MPI_COMM_WORLD, &status);
/* Send edge data to our left neighbour, receive from right.
We are sending the edge (locV[1][1]..locV[locNmyo][1]),
and receiving into edge (locV[1][locTmyo+1]..locV[locNmyo][locTmyo+1]) */
MPI_Sendrecv(&(locV[1][1]), 1, leftright, left, lefttag,
&(locV[1][locTmyo+1]), 1, leftright, right, lefttag,
MPI_COMM_WORLD, &status);
/* Send edge data to our up neighbour, receive from down.
We are sending the edge (locV[1][1]..locV[1][locTmyo]),
and receiving into edge (locV[locNmyo+1][1]..locV[locNmyo+1][locTmyo]) */
MPI_Sendrecv(&(locV[1][1]), 1, updown, up, uptag,
&(locV[locNmyo+1][1]), 1, updown, down, uptag,
MPI_COMM_WORLD, &status);
/* Send edge data to our down neighbour, receive from up.
We are sending the edge (locV[locNmyo][1]..locV[locNmyo][locTmyo]),
and receiving into edge (locV[0][1]..locV[0][locTmyo]) */
MPI_Sendrecv(&(locV[locNmyo][1]),1, updown, down, downtag,
&(locV[0][1]), 1, updown, up, downtag,
MPI_COMM_WORLD, &status);
/* Release the resources associated with the Type_create() calls. */
MPI_Type_free(&updown);
MPI_Type_free(&leftright);
}
void printArrays(char **locV, const int locNmyo, const int locTmyo,
const int size, const int rank) {
/* all these barriers are a terrible idea, but it's just
for controlling output to the screen as a demo. You'd
really do something smarter here... */
for (int task=0; task<size; task++) {
if (rank == task) {
printf("\nTask %d's local array:\n", rank);
for (int i=0; i<locNmyo+2; i++) {
putc('[', stdout);
for (int y=0; y<locTmyo+2; y++) {
putc(locV[i][y], stdout);
}
printf("]\n");
}
}
fflush(stdout);
MPI_Barrier(MPI_COMM_WORLD);
}
}
int main(int argc, char **argv) {
int ierr, size, rank;
char **locV;
const int Nmyo=12; /* horizontal */
const int Tmyo=12; /* vertical */
const int ncols=2; /* n procs in horizontal direction */
int nrows;
int myrow, mycol;
int locNmyo, locTmyo;
ierr = MPI_Init(&argc, &argv);
ierr|= MPI_Comm_size(MPI_COMM_WORLD, &size);
ierr|= MPI_Comm_rank(MPI_COMM_WORLD, &rank);
nrows = size/ncols;
if (nrows*ncols != size) {
fprintf(stderr,"Size %d does not divide number of columns %d!\n",
size, ncols);
MPI_Abort(MPI_COMM_WORLD,-1);
}
/* where are we? */
mycol = rank % ncols;
myrow = rank / ncols;
/* figure out how many Tmyo we have */
locTmyo = (Tmyo / ncols);
/* in case it doesn't divide evenly... */
if (mycol == ncols-1) locTmyo = Tmyo - (ncols-1)*locTmyo;
/* figure out how many Tmyo we have */
locNmyo = (Nmyo / nrows);
/* in case it doesn't divide evenly... */
if (myrow == nrows-1) locNmyo = Nmyo - (ncols-1)*locNmyo;
/* allocate our local array, with space for edge data */
locV = alloc_2d_char(locNmyo+2, locTmyo+2);
/* fill in our local data - first spaces everywhere */
for (int i=0; i<locNmyo+2; i++)
for (int y=0; y<locTmyo+2; y++)
locV[i][y] = ' ';
/* then the inner regions have our rank # */
for (int i=1; i<locNmyo+1; i++)
for (int y=1; y<locTmyo+1; y++)
locV[i][y] = '0' + rank;
/* The "before" picture: */
if (rank==0) printf("###BEFORE###\n");
printArrays(locV, locNmyo, locTmyo, size, rank);
/* Now do edge filling. Ignore corners for now;
the right way to do that depends on your algorithm */
edgeDataFill(locV, locNmyo, locTmyo, ncols, myrow, mycol, size, rank);
/* The "after" picture: */
if (rank==0) printf("###AFTER###\n");
printArrays(locV, locNmyo, locTmyo, size, rank);
MPI_Finalize();
}
Вышеприведенную программу можно еще больше упростить, используя MPI_Cart_create, чтобы создать многомерный домен и автоматически рассчитать для вас соседей, но я хотел показать вам логику, чтобы вы видели, что происходит.
Кроме того, если вы можете посоветоваться с кем-то, кто делал это долгое время:
Каждый раз, когда у вас есть строка за строкой повторяющегося кода: примерно 60 (!!) строк этого:
Vmax =V[i][y]-Vold; updateMaxStateChange(Vmax / dt);
mmax=m[i][y]-mold; updateMaxStateChange(mmax / dt);
hmax=h[i][y]-hold; updateMaxStateChange(hmax / dt);
jmax=j[i][y]-jold; updateMaxStateChange(jmax / dt);
mLmax=mL[i][y]-mLold; updateMaxStateChange(mLmax / dt);
hLmax=hL[i][y]-hLold; updateMaxStateChange(hLmax / dt);
hLBmax=hLB[i][y]-hLBold; updateMaxStateChange(hLBmax / dt);
hLSmax=hLS[i][y]-hLSold; updateMaxStateChange(hLSmax / dt);
amax=a[i][y]-aold; updateMaxStateChange(amax / dt);
i1fmax=i1f[i][y]-i1fold; updateMaxStateChange(i1fmax / dt);
i1smax=i1s[i][y]-i1sold; updateMaxStateChange(i1smax / dt);
Xrmax=Xr[i][y]-Xrold; updateMaxStateChange(Xrmax / dt);
i2max=i2[i][y]-i2old; updateMaxStateChange(i2max / dt);
это признак того, что вы не используете правильные структуры данных. Здесь вы почти наверняка захотите иметь трехмерный массив переменных состояния, при этом (вероятно) 3-й индекс - это видовая или локальная переменная состояния или как вы хотите назвать i2, i1f, i1s и т. Д. Затем все эти строки можно заменить с циклом, и добавление новой локальной переменной состояния становится намного проще.
Точно так же, когда все ваше состояние будет определено как глобальные переменные, это сделает вашу жизнь намного сложнее, когда дело доходит до обновления и обслуживания кода. Опять же, это, вероятно, отчасти связано с наличием вещей в миллионах независимых переменных состояния вместо того, чтобы иметь структуры или массивы более высокой размерности, объединяющие все соответствующие данные.