c - 在我的possion方程代码中集体中止

Question

/***every function is working correct but after only first iteration is giving collective abort anyone can tell what is or coulde be the reason***/
#include<stdio.h>
#include<stdlib.h>
#include"mpi.h"

const double tolerance = 0.00001;
const int maxit = 10000;

void MPE_decomp1d(int n, int size, int id, int *s, int *e)
{
    /*****calculating start and end row for every process*****/
    *s = (n/size)*id + ((n%size)>0)*(id>(n%size)?n%size:id);
        *e = *s + (n/size)+((n%size)>id);
}

void onedinit(double **a, double **b, double **f, const int nx, const int s, const int e)
{
    int i, j;
    int ls, le;
    ls = s - (s!=0); 
    le = e + (e!=nx);
    /***setting all the intial values to zero****/
    for (i = ls; i < le; i++)
    {
        for (j = 0; j < nx; j++)
        {
            a[i][j] = b[i][j] = f[i][j] = 0.0;
        }
    }
    //***************************Defining Boundary Condition***********************************//
    /***setting left boundary to 1***/
    for (i = ls; i < le; i++) a[i][0] = b[i][0] = 1;
    /***setting value f(0, i) is 2***/
    if (s==0) for (i = 0; i < nx; i++) a[0][i] = b[0][i] = 2.0;
}   

void exchng1(double **a, const int nx, const int s, const int e, MPI_Comm comm1d, int nbrbottom, int nbrtop)
{
    int rank, coord;
    MPI_Status status;
    MPI_Comm_rank(comm1d, &rank);
    MPI_Cart_coords(comm1d, rank, 1, &coord);
    /*****************if process id is odd then first send and if even then first recive to avoid deadlock**********/
    if (coord&1)
    {
        if (nbrbottom != -1) MPI_Send(a[e-s], nx, MPI_DOUBLE, nbrbottom, 0, comm1d);  
        if (nbrtop != -1)    MPI_Recv(a[0], nx, MPI_DOUBLE, nbrtop, 1, comm1d, &status);
        if (nbrtop != -1)    MPI_Send(a[1], nx, MPI_DOUBLE, nbrtop, 0, comm1d);
        if (nbrbottom != -1) MPI_Recv(a[e-s+1], nx, MPI_DOUBLE, nbrbottom, 1, comm1d, &status);
    }
    else
    {
        if (nbrtop != -1)    MPI_Recv(a[0], nx, MPI_DOUBLE, nbrtop, 0, comm1d, &status);
        if (nbrbottom != -1) MPI_Send(a[e-s-(s==0)], nx, MPI_DOUBLE, nbrbottom, 1, comm1d);
        if (nbrbottom != -1) MPI_Recv(a[e-s+(s!=0)], nx, MPI_DOUBLE, nbrbottom, 0, comm1d, &status);
        if (nbrtop != -1)    MPI_Send(a[1], nx, MPI_DOUBLE, nbrtop, 1, comm1d);
    }
}

void sweep1d(double **a, double **f, int nx, const int s, const int e, double **b)
{
    int i, j;
    int rows;
    rows = e - s - (s==0) - (e==0);
    nx -= 1;
    double h = 1.0 / (double)nx;
    for (i = 1; i <= rows; i++) for (j = 1; j < nx; j++)
        b[i][j] = 0.25 * (a[i-1][j] + a[i][j+1] + a[i][j-1] + a[i+1][j]) - h*h*f[i][j];
    return;
}

double diff(double **a, double **b, const int nx, int s, int e)
{
    double sum = 0.0;
    int i, j;
    int st, ed;
    st = (s!=0);
    ed = e-s+(s!=0);
    for (i = st; i < ed; i++) for (j = 0; j < nx; j++)
        sum += (a[i][j] - b[i][j])*(a[i][j] - b[i][j]);
    return sum;
}
int main(int argc, char *argv[])
{
    int nx, ny;
    int myid, root, numprocs, period=0;
    int nbrbottom, nbrtop, s, e, it;

    double diffnorm, dwork;
    double t1, t2;
    double **a, **b, **f;

    root = 0;
    MPI_Comm comm1d;
    MPI_Init(&argc, &argv);;
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    if(!myid)
    {
        /******for this piece of code nx and ny are assumed to be same please*******/
        printf("Enter the number of cells in X & Y direction\n");
        scanf("%d %d", &nx, &ny);
        nx += 1;
        ny += 1;
        ny = nx; //forced to follow our assumption;
    }
    MPI_Bcast(&nx, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&ny, 1, MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Cart_create(MPI_COMM_WORLD, 1, &numprocs, &period, 1, &comm1d);

    MPI_Comm_rank(comm1d, &myid);
    MPI_Cart_shift(comm1d, 0, 1, &nbrtop, &nbrbottom);

    MPE_decomp1d(ny, numprocs, myid, &s, &e);
    int ls, le, rows;
    int i, j;
    ls = s - (s!=0); 
    le = e + (e!=nx);
    rows = le - ls;

    a = (double**)malloc(rows*sizeof(double*));
    b = (double**)malloc(rows*sizeof(double*));
    f = (double**)malloc(rows*sizeof(double*));
    for (i = ls; i < le; i++)
    {
        a[i] = (double*)malloc(nx*sizeof(double));
        b[i] = (double*)malloc(nx*sizeof(double));
        f[i] = (double*)malloc(nx*sizeof(double));
    }
    onedinit(a, b, f, nx, s, e);
    diffnorm = 0.0;
    it = 0;
    do
    {
//      printf("%danshu\n", myid);

        exchng1(a, nx, s, e, comm1d, nbrbottom, nbrtop);
        sweep1d(a, f, nx, s, e, b);
        exchng1(b, nx, s, e, comm1d, nbrbottom, nbrtop);
        sweep1d(b, f, nx, s, e, a);
        dwork = diff(a, b, nx, s, e);

        /************printing matrix a after every iteration******/
        for (i = 0; i < rows; i++) 
        {
            for (j = 0; j < nx; j++) printf("%lf ", a[i][j]);
            printf("\n");
        }

        MPI_Barrier(comm1d);

        //printf("%lfhehe\n", dwork);
        MPI_Allreduce(&dwork, &diffnorm, 1, MPI_DOUBLE, MPI_SUM, comm1d);
        //printf("%dhere\n", myid);
    } 
    while (++it < maxit && diffnorm > tolerance);

    MPI_Finalize();
    return 0;
}

Answer 1

因此，仅在 SO 中转储 130 行代码并询问为什么它不起作用可能不是获得良好答案的最佳方式 - 特别是当您编写的唯一实际句子是“每个函数都在工作”时......如果那样的话，你就不会有问题了。您需要将范围缩小到更具体的案例，并提出更具体的问题。

在这种特殊情况下，我在过去的教学中看到过很多这样的代码，因此可以了解一些正在发生的事情。

首先，你不能做这样的事情:

ls = s - (s!=0); 
le = e + (e!=nx);
rows = le - ls;

a = (double**)malloc(rows*sizeof(double*));
/*...*/
for (i = ls; i < le; i++)
{
    a[i] = (double*)malloc(nx*sizeof(double));
    /*...*/
}

如果您将 100 行分成 4 个处理器，并且您是(比如说)MPI 任务 2，那么您的 s 是 50，e 是 75，并且所以 ls 会是 49 而 le 会是 76，所以你试图访问 a[49..76] 即使你只分配了大小为 27 的 a!该特定错误出现在整个代码中，需要修复。您想要访问 a[0..rows-1]。

顺便说一下，我什至没有检查过 MPE_decomp1d 是否真的做了正确的事情。我们都经历过这样一个阶段，我们认为在 C 中通过使用逻辑表达式乘以三元运算符等将内容放在一行中很可爱，但严重的是，当其他人必须修复它时，它会使您的代码变得不必要的乏味 - 2个月后无论是SOers还是你自己。

在 exchng1 中，您正在做不必要的工作。您不需要检查 nbrbottom 或 nbrtop 是否有效；如果不是，MPI_Cart_shift 返回 MPI_PROC_NULL，发送或接收是空操作。所以从这些队列发送/接收是无害的，这是一个很好的设计决策，因为它避免了逻辑中的许多极端情况。

同样，为避免死锁，您可以使用 MPI_Sendrecv 而不是单独的 Send 和 Recv。再加上上面的意思，而不是这个:

if (coord&1)
{
    if (nbrbottom != -1) MPI_Send(a[e-s], nx, MPI_DOUBLE, nbrbottom, 0, comm1d);  
    if (nbrtop != -1)    MPI_Recv(a[0], nx, MPI_DOUBLE, nbrtop, 1, comm1d, &status);
    if (nbrtop != -1)    MPI_Send(a[1], nx, MPI_DOUBLE, nbrtop, 0, comm1d);
    if (nbrbottom != -1) MPI_Recv(a[e-s+1], nx, MPI_DOUBLE, nbrbottom, 1, comm1d, &status);
}
else
{
    if (nbrtop != -1)    MPI_Recv(a[0], nx, MPI_DOUBLE, nbrtop, 0, comm1d, &status);
    if (nbrbottom != -1) MPI_Send(a[e-s-(s==0)], nx, MPI_DOUBLE, nbrbottom, 1, comm1d);
    if (nbrbottom != -1) MPI_Recv(a[e-s+(s!=0)], nx, MPI_DOUBLE, nbrbottom, 0, comm1d, &status);
    if (nbrtop != -1)    MPI_Send(a[1], nx, MPI_DOUBLE, nbrtop, 1, comm1d);
}

你可以这样做:

MPI_Sendrecv(a[e-s], nx, MPI_DOUBLE, nbrbottom, 0, a[0], nx, MPI_DOUBLE, nbrtop, 0, comm1d, &status);
MPI_Sendrecv(a[1], nx, MPI_DOUBLE, nbrtop, 1, a[e-s+1], nx, MPI_DOUBLE, nbrbottom, 1, comm1d, &status);

-- 更简单，对吧？

不过交换还是有一些问题；接收到 a[e-s+1] 是不对的，尽管正如我所提到的，我懒得解密 MPE_decomp1d 来找出原因。假设您希望接收到 a[rows-1]。

最后，MPI_Barrier() 很慢而且完全没有必要；在 guardcell 交换中有足够的同步(更不用说 Allreduce 了)你不需要它。

完成所有这些更改后，代码运行时不会出现内存访问问题；你必须检查它是否给出了正确的答案。

#include<stdio.h>
#include<stdlib.h>
#include"mpi.h"

const double tolerance = 0.00001;
const int maxit = 10000;

void MPE_decomp1d(int n, int size, int id, int *rows)
{
    int s, e;
    s = (n/size)*id + ((n%size)>0)*(id>(n%size)?n%size:id);
    e = s + (n/size)+((n%size)>id);
    *rows = e - s - (s==0) - (e==0);
}

void onedinit(double **a, double **b, double **f, const int nx, const int rows, const int id, const int nprocs)
{
    int i, j;

    for (i = 0; i < rows; i++)
    {
        for (j = 0; j < nx; j++)
        {
            a[i][j] = b[i][j] = f[i][j] = 0.0;
        }
    }
    for (i = 0; i < rows; i++) a[i][0] = b[i][0] = 1;

    if (id == 0)
        for (i = 0; i < nx; i++) a[0][i] = b[0][i] = 2.0;
}
void exchng1(double **a, const int nx, const int rows, MPI_Comm comm1d, int nbrbottom, int nbrtop)
{
    int rank, coord;
    MPI_Status status;
    MPI_Comm_rank(comm1d, &rank);
    MPI_Cart_coords(comm1d, rank, 1, &coord);

    /* send data downwards */
    MPI_Sendrecv(a[rows-2], nx, MPI_DOUBLE, nbrbottom, 0, a[0], nx, MPI_DOUBLE, nbrtop, 0, comm1d, &status);
    /* send data upwards */
    MPI_Sendrecv(a[1], nx, MPI_DOUBLE, nbrtop, 1, a[rows-1], nx, MPI_DOUBLE, nbrbottom, 1, comm1d, &status);
}

void sweep1d(double **a, double **f, const int nx, const int rows, double **b)
{
    int i, j;
    double h = 1.0 / (double)nx;
    for (i = 1; i < rows-1; i++) for (j = 1; j < nx-1; j++)
        b[i][j] =  
            0.25 * ( a[i-1][j] + a[i][j+1] + a[i][j-1] + a[i+1][j]) - h*h*f[i][j];
    return;
}

double diff(double **a, double **b, const int nx, const int rows)
{   
    double sum = 0.0;
    int i, j;
    for (i = 0; i < rows; i++) for (j = 0; j < nx; j++)
        sum += (a[i][j] - b[i][j])*(a[i][j] - b[i][j]);
    return sum;
}
int main(int argc, char *argv[])
{   
    int nx, ny;
    int myid, root, numprocs, period=0;
    int nbrbottom, nbrtop, it;

    double diffnorm, dwork;

    double **a, **b, **f;

    root = 0;
    MPI_Comm comm1d;
    MPI_Init(&argc, &argv);;
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

    if(!myid)
    {   
        /******for this piece of code nx and ny are assumed to be same please*******/
        printf("Enter the number of cells in X & Y direction\n");
        scanf("%d %d", &nx, &ny);
        nx += 1;
        ny += 1; 
        ny = nx; //forced to follow our assumption;
    }
    MPI_Bcast(&nx, 1, MPI_INT, 0, MPI_COMM_WORLD);
    MPI_Bcast(&ny, 1, MPI_INT, 0, MPI_COMM_WORLD);

    MPI_Cart_create(MPI_COMM_WORLD, 1, &numprocs, &period, 1, &comm1d);

    MPI_Comm_rank(comm1d, &myid);
    MPI_Cart_shift(comm1d, 0, 1, &nbrtop, &nbrbottom);

    int rows;
    MPE_decomp1d(ny, numprocs, myid, &rows);
    int i, j;

    a = (double**)malloc(rows*sizeof(double*));
    b = (double**)malloc(rows*sizeof(double*));
    f = (double**)malloc(rows*sizeof(double*));
    for (i = 0; i < rows; i++)
    {
        a[i] = (double*)malloc(nx*sizeof(double));
        b[i] = (double*)malloc(nx*sizeof(double));
        f[i] = (double*)malloc(nx*sizeof(double));
    }
    onedinit(a, b, f, nx, rows, myid, numprocs);
    diffnorm = 0.0;
    it = 0;
    do
    {
        exchng1(a, nx, rows, comm1d, nbrbottom, nbrtop);
        sweep1d(a, f, nx, rows, b);
        exchng1(b, nx, rows, comm1d, nbrbottom, nbrtop);
        sweep1d(b, f, nx, rows, a);
        dwork = diff(a, b, nx, rows);

        /************printing matrix a after every iteration******/
        for (i = 0; i < rows; i++)
        {
            for (j = 0; j < nx; j++) printf("%lf ", a[i][j]);
            printf("\n");
        }

        //printf("%lfhehe\n", dwork);
        MPI_Allreduce(&dwork, &diffnorm, 1, MPI_DOUBLE, MPI_SUM, comm1d);
        //printf("%dhere\n", myid);
    }
    while (++it < maxit && diffnorm > tolerance);

    MPI_Finalize();
    return 0;
}