mpi并行求解二维泊松方程，编译连接没有问题，结果出错，希望指教。

Question

#include "stdafx.h"
#define MPICH_SKIP_MPICXX
#include

#include

#include

#include

#include

#pragma comment (lib, "mpi.lib")

#include "mpi.h"

#include

double max(double a,double b)
{
return a>=b?a:b;

}

double dw=2.0,dh=3.0;//求解区域的大小
int im=30,jm=60;//全局网格规模
int npx=1,npy=1;//XY方向进程数目
const iml=im/npx,jml=jm/npy;//各进程局部网格规模
double u[31][61];//定义网格节点近似解
double us[31][61];//定义网格节点精确解
double u0[31][61];//定义Jacobi迭代辅助变量
double f[31][61],solution[31][61],rhs[31][61];//f(x,y)在各节点的值
int nproc;
int main(int argc, char* argv[])

{

int myrank,myleft,myright,myup,mydown;//各进程的进程号，四个相邻进程的进程号
int mypx,mypy;//各进程的进程号沿X,Y方向的坐标
double xst,yst;//各进程拥有的子区域沿X,Y方向的起始坐标
double hx,hy;//X,Y方向离散网格步长
double hx2,hy2,hxy2,rhxy;
int ist,iend,jst,jend;//各进程沿X，Y方向内部网格节点的起始和终止坐标
int htype,vtype;//MPI用户自定义数据类型，表示各进程沿X,Y方向与相邻进程交换数据单元
double t0,t1;
for(int i=0;i<=iml;i++)
{
for(int j=0;j<=jml;j++)
{
solution[i][j]=i*i+i*j;//解析解
rhs[i][j]=-4.0;//泊松方程源项
u[i][j]=0;
}
}
int namelen;

char processor_name[MPI_MAX_PROCESSOR_NAME];

MPI_Init(&argc, &argv);   
//用MPI_Comm_size 获得进程个数
int size;
size=npx*npy;//进程数必须等于NPX*NPY
MPI_Comm_size(MPI_COMM_WORLD,&nproc);
struct MPI_Status *status = (   struct MPI_Status *)malloc(size * sizeof(int));
if(nproc!=size||(im%npx!=0)||(jm%npy!=0))
{
  printf("error!exit out");
  MPI_Finalize(); 
  MPI_Get_processor_name(processor_name, &namelen);    
  MPI_Finalize();
} 
//按自然顺序确定各进程自身及其相邻四个进程的进程号
MPI_Comm_rank(MPI_COMM_WORLD, &myrank); 
myleft=myrank-1;
if(myrank%npx==0){myleft=MPI_PROC_NULL;
myright=myrank+1;}
if(myright%npx==0){myright=MPI_PROC_NULL;
myup=myrank+npx;}
if(myup>=size){myup=MPI_PROC_NULL;
mydown=myrank-npx;}
if(mydown<0){mydown=MPI_PROC_NULL;
}
mypy=myrank/npx;
mypx=myrank-mypy*npx;
hx=dw/im;
hx2=hx*hx;
hy=dh/jm;
hy2=hy*hy;
hxy2=hx2*hy2;
rhxy=0.5/(hx2+hy2);
double dx=hx2*rhxy;
double dy=hy2*rhxy;
double dd=rhxy*hxy2;
xst=mypx*dw/npx;
yst=mypy*dh/npy;
ist=1;
iend=iml;
if(mypx==(npx-1)){iend=iend-1;}//最右边区域X方向少一个点
jst=1;
jend=jml;
if(mypy==(npy-1)){jend=jend-1;}//最上边区域Y方向少一个点
MPI_Type_contiguous(iend-ist+1,MPI_DOUBLE,&htype);
MPI_Type_commit(&htype);
MPI_Type_vector(jend-jst+1,1,iml+2,MPI_DOUBLE,&vtype);
MPI_Type_commit(&vtype);
for(int j=jst-1;j<=jend+1;j++)
{
   for(i=ist-1;i<=iend+1;i++)
   {
      int xx=(int)((i+mypx*iml)*hx);
      int yy=(int)((j+mypy*jml)*hy);
      if(i>=ist&&i<=iend&&j>=jst&&j<=jend){
      u[i][j]=0;
      us[i][j]=solution[xx][yy];
      f[i][j]=dd*rhs[xx][yy];}
      if((i==ist-1&&mypx==0)||(j==jst-1&&mypy==0)||(i==iend+1&&mypx==npx-1)||(j==jend+1&&mypy==npy-1)){
      u[i][j]=solution[xx][yy];}
   }
}
//Jacobi迭代
int niter=0;
t0=MPI_Wtime();

tt:niter=niter+1;
//交换辅助网格上的近似解
MPI_Send(&u[1][1],1,vtype,myleft,niter+100,MPI_COMM_WORLD);//发送左边界
//int MPI_Send(void buf,int count,MPI_Datatype datatype,int dest,int tag,MPI_Comm comm);
MPI_Send(&u[iend][1],1,vtype,myright,niter+100,MPI_COMM_WORLD);//发送右边界
MPI_Send(&u[1][1],1,htype,mydown,niter+100,MPI_COMM_WORLD);//发送下边界
MPI_Send(&u[1][jend],1,htype,myup,niter+100,MPI_COMM_WORLD);//发送上边界
MPI_Recv(&u[iend][1],1,vtype,myright,niter+100,MPI_COMM_WORLD,status);//接收右边界
MPI_Recv(&u[0][1],1,vtype,myleft,niter+100,MPI_COMM_WORLD,status);//接收左边界
MPI_Recv(&u[1][jend],1,htype,myup,niter+100,MPI_COMM_WORLD,status);//接收上边界
MPI_Recv(&u[1][0],1,htype,mydown,niter+100,MPI_COMM_WORLD,status);//接收下边界
for(j=jst;j<=jend;j++)
{
for(i=ist;i<=iend;i++)
{
u0[i][j]=f[i][j]+dx(u[i][j-1]+u[i][j+1])+dy*(u[i-1][j]+u[i+1][j]);
}
}
//计算误差
double err=0.0;
for(j=jst;j<=jend;j++)
{
for(i=ist;i<=iend;i++)
{
u[i][j]=u0[i][j];
err=max(err,fabs(u[i][j]-us[i][j]));
}
}
double err0=err;
MPI_Reduce(&err0,&err,1,MPI_DOUBLE,MPI_MAX,0,MPI_COMM_WORLD);
if(myrank==0&&niter%100==0){
printf("niter=%f,err=%f",niter,err);}
if(err>=5){goto tt;}
t1=MPI_Wtime();
if(myrank==0){
printf("successfully converged after %f iterations,error=%f,wtime=%f",niter,err,t1-t0);}
MPI_Finalize();
free(status);
return 0;

}