/*******************************************************************************
Erik Erhardt
Math 471
Project 3
heat.c

on frontend:
  make heat
  mpirun -np 25 ./heat 1 400 5 5
  mpirun -np 16 ./heat 2 400 4 4

on backend (many nodes):
  qsub -I -l nodes=2,walltime=900
  mpirun -np 4 -machinefile $PBS_NODEFILE ./heat 2 100 2 2
  qsub -I -l nodes=4,walltime=900
  mpirun -np 16 -machinefile $PBS_NODEFILE ./heat 2 200 4 4
  qsub -I -l nodes=7,walltime=900
  mpirun -np 25 -machinefile $PBS_NODEFILE ./heat 2 200 5 5

  qsub -I -l nodes=1,walltime=900
  mpirun -np 1 ./heat 1 100 1 1

%%%%%%%%%%%%%%%% part 1 laplace
  mpirun -np 1  ./heat 1 100 1 1
  mpirun -np 1  ./heat 1 200 1 1
  mpirun -np 1  ./heat 1 400 1 1
  mpirun -np 1  ./heat 1 800 1 1
  mpirun -np 25 ./heat 1 100 5 5
  mpirun -np 25 ./heat 1 200 5 5
  mpirun -np 25 ./heat 1 400 5 5
  mpirun -np 25 ./heat 1 800 5 5
  mpirun -np 25 ./heat 1 100 1 25
  mpirun -np 25 ./heat 1 200 1 25
  mpirun -np 25 ./heat 1 400 1 25
  mpirun -np 25 ./heat 1 800 1 25

         1x1           5x5           1x25
100x100  0.0015910596  0.0015910596  0.0015910596
200x200  0.0004018149  0.0004018149  0.0004018149
400x400  0.0001009606  0.0001009606  0.0001009606
800x800  0.0000253035  0.0000253035  0.0000253035


%%%%%%%%%%%%%%%% part 2 scalability

  qsub -I -l nodes=7,walltime=900
  mpirun -np  1 -machinefile $PBS_NODEFILE ./heat 3 200 1 1
  mpirun -np  4 -machinefile $PBS_NODEFILE ./heat 3 200 2 2
  mpirun -np 16 -machinefile $PBS_NODEFILE ./heat 3 200 4 4
  mpirun -np 25 -machinefile $PBS_NODEFILE ./heat 3 200 5 5
  mpirun -np  1 -machinefile $PBS_NODEFILE ./heat 3 200 1  1
  mpirun -np  4 -machinefile $PBS_NODEFILE ./heat 3 200 1  4
  mpirun -np 16 -machinefile $PBS_NODEFILE ./heat 3 200 1 16
  mpirun -np 25 -machinefile $PBS_NODEFILE ./heat 3 200 1 25


%%%%%%%%%%%%%%%% part 2 contour
  qsub -I -l nodes=7,walltime=900
  mpirun -np 25 -machinefile $PBS_NODEFILE ./heat 2 800 5 5

  qsub -I -l nodes=1,walltime=900
  mpirun -np 1 ./heat 2 100 1 1








*******************************************************************************/

#include <mpi.h>
#include <math.h>
#include <stdio.h>

#define MAX_DIM 800

int rank, rank_x, rank_y, rank_down, rank_up, rank_left, rank_right;
int nprocs, nprocsx, nprocsy;
int nxg, nyg, nx, ny, n;
double length_x, length_y, delta_x, delta_y;
double PI, h, t0, tfinal, t, time_per_step;
int junk, error_flag=0, sw_method;
double U[MAX_DIM+2][MAX_DIM+2], Ue[MAX_DIM+2][MAX_DIM+2], lu[MAX_DIM+2][MAX_DIM+2];  // this is our data matrix


/******************************************************************************/
int main(int argc, char *argv[]){
   int i, j, it;
   double time1, time2, timetotal, timesofar, timetogo;          // time the run

   MPI_Init(&argc,&argv);                   // initialize MPI
   MPI_Comm_rank(MPI_COMM_WORLD, &rank);    // my rank number
   MPI_Comm_size(MPI_COMM_WORLD, &nprocs);  // number of processes

   time1 = MPI_Wtime();                     // first time t1

   junk = init(argc, &*argv);  /*******/                // initialize values

   if(sw_method==2){
      for (it=0; it<=n+1; it++){
         if(it % ((n-1)/100) == 0)
         { // output grid to a file
            double x_out[nx];
            double y_out[ny];
            double u_out[nx][ny];

            for (i=0; i<nx; ++i)   x_out[i]=( i + (rank_x*nx))/((double)(nxg-1));
            for (j=0; j<ny; ++j)   y_out[j]=( j + ((nprocsy-rank_y-1)*ny))/((double)(nyg-1));

            for (i=0; i<nx; i++) {
               for (j=0; j<ny; j++) {
                  u_out[i][j]=U[i+1][j+1];
               }
            }
            junk = output(nx,ny,x_out,y_out,u_out,t);
         }
         if((rank==0) && (it%10000==0)) {
            time2 = MPI_Wtime();                         // second time t2
            timesofar = time2 - time1;                   // total time is difference
            timetogo = (timesofar/it)*(n-it);
            printf("it=%d/%d, t=%f (%2.1f\%) (realtime so far=%f, to go=%f)\n", it, n, t, (double)100*it/n, timesofar, timetogo);
         }
         junk = euler();
         t=t+h;
      }
   }

   if(sw_method==3){
      time1 = MPI_Wtime();                     // reset time for method 3 first time t1
      for (it=0; it<=n+1; it++){
         junk = euler();
         t=t+h;
      }
   }

   time2 = MPI_Wtime();                               // second time t2
   timetotal = time2 - time1;                         // total time is difference
   time_per_step = timetotal/(n+2);
   if (rank==0) printf("time=%f, time_per_step=%f\n", timetotal, time_per_step);
   //if (rank==0) printf("time=%f\n", timetotal);

   MPI_Finalize();
}


/******************************************************************************/
/* init *****************************************/
int init(int argc, char **argv){
   int i, j;
   int argc_counter, string_counter;
   int result;
   char argv1[81];

   if(argc<5){
       if(rank==0) printf("argc %d\n", argc);
       if(rank==0) printf("Command line args: sw_method(1,2) nxg nprocsx nprocsy\n");
       if(rank==0) printf("EXITING\n");
       exit(0);
   }

   // sw_method=1 is part 1 of the assignment testing Laplacian,
   //           2 is part 2 for propagation of Marshak wave.
   strcpy(argv1,argv[1]); sw_method = string2int(argv1);
   strcpy(argv1,argv[2]); nxg       = string2int(argv1);
   strcpy(argv1,argv[3]); nprocsx   = string2int(argv1);
   strcpy(argv1,argv[4]); nprocsy   = string2int(argv1);

   if(rank==0) printf("Params: sw_method=%d, nxg=%d, nprocsx=%d, nprocsy=%d\n", sw_method, nxg, nprocsx, nprocsy);

   PI=4*atan( (double)1 );                  // known value of pi
   length_x=1; length_y=length_x;
   nyg=nxg;
   delta_x=length_x/(nxg+1); delta_y=delta_x;

   if(nxg%nprocsx==0){nx=nxg/nprocsx;}else{if(rank==0){printf("ERROR: nxg %% nprocsx!=0 (%d %d)\n", nxg, nprocsx);}  error_flag++;}
   if(nyg%nprocsy==0){ny=nyg/nprocsy;}else{if(rank==0){printf("ERROR: nyg %% nprocsy!=0 (%d %d)\n", nyg, nprocsy);}  error_flag++;}

   if(nprocs != nprocsx*nprocsy){if(rank==0){printf("ERROR: nprocs != nprocsx*nprocsy (%d != %d * %d), EXIT\n", nprocs, nprocsx, nprocsy);} error_flag++;}

   if(sw_method == 1){
      if(rank==0) printf("Method 1, testing Laplacian.\n");
      n=1;
   };

   if(sw_method==2 | sw_method==3){
      if(rank==0) printf("Method 2, Marshak wave.\n");
      h=0.0000001; t0=0.0; tfinal=0.5; t=t0;
      h=(pow(delta_x,2)/8)/4;
      //h=0.0001; t0=0.0; tfinal=0.01; t=t0;
      n = ceil((tfinal-t0)/h)+1;       // number of timesteps
      if(sw_method==3){ n=100;}
      if(rank==0) printf("Params: h=%f, t0=%f, tfinal=%f, n=%d\n", h, t0, tfinal, n);
   }

   if(h >= pow(delta_x,2)/8){if(rank==0){printf("ERROR: h >= pow(delta_x,2)/8 (%e >= %e ), EXIT\n", h, pow(delta_x,2)/8);} error_flag++;}

   if(error_flag>0){if(rank==0) printf("Errors = %d\n EXITING\n", error_flag); exit(0);}

   junk = init_local_ranks();
   junk = init_block_data();

   return 0;
}

int init_local_ranks(){
   // if(rank==0) printf("nprocs x y %d %d %d\n", nprocs, nprocsx, nprocsy);

   // x, y rank location of current process
   rank_x = rank % nprocsx;   // rank in x direction (mod)
   rank_y = rank / nprocsx;   // rank in y direction (int division)

   /* Ranks (blocks) go in this order
      6 7 8
      3 4 5
      0 1 2  */

   // ranks of neighboring processors
   rank_down  = (rank_y-1)*nprocsx + rank_x  ;    if(rank_y  ==  0       ) rank_down  = -1;
   rank_up    = (rank_y+1)*nprocsx + rank_x  ;    if(rank_y  == nprocsy-1) rank_up    = -1;
   rank_left  =  rank_y   *nprocsx + rank_x-1;    if(rank_x  ==  0       ) rank_left  = -1;
   rank_right =  rank_y   *nprocsx + rank_x+1;    if(rank_x  == nprocsx-1) rank_right = -1;

   //printf("Ranks: me, x, y, up, down, left, right: %d %d %d %d %d %d %d \n", rank, rank_x, rank_y, rank_up, rank_down, rank_left, rank_right);
   return 0;
}

int init_block_data(){
   int i, j;
   double b, initx[nx+2], inity[nx+2], pi2;

   /* Data goes in this order per block (border is ghost cells):
      0 1 2 3 4 ... nx+1
      1
      2    ...
      ...
      ny+1             */

   for (i=0; i<nx+2; i++){for (j=0; j<ny+2; j++){U[i][j]=0;}}  // init at 0

   if(sw_method==1){
      pi2 = pow(PI,2);
      //printf("length_x,length_y,nxg,nyg,nx,ny,rank_x,rank_y %f %f %d %d %d %d %d %d\n",length_x,length_y,nxg,nyg,nx,ny,rank_x,rank_y);
      for(i=0; i<nx+2; i++){initx[i] = (length_x/(nxg+1))*(nx*rank_x+i);}
      for(j=0; j<ny+2; j++){inity[j] = (length_y/(nyg+1))*(ny*(nprocsy-1-rank_y)+j);}
      for(j=0; j<ny+2; j++){
         for(i=0; i<nx+2; i++){
             U[i][j] =        sin(PI*initx[i])*sin(PI*inity[j]);  // test function
            Ue[i][j] = -2*pi2*sin(PI*initx[i])*sin(PI*inity[j]);  // test function (exact value)
            //printf("rank=%d i=%3d j=%3d %10.5f %10.5f %10.5f %10.5f\n",rank,i,j,initx[i],inity[j],U[i][j],Ue[i][j]);
         }
      }
      junk = laplacian();
      junk = max_error();
   }

   if(sw_method==2){
      // init ghost cell boundary conditions
      if(rank_down==-1){
         for(i=0; i<nx+2; i++){
            initx[i] = (length_x/(nxg+1))*(nx*rank_x+i);
            inity[i] = pow(sin(PI*initx[i]),2);
            // printf("rank = %d, initx %d %f %f\n", rank, i, initx[i], inity[i]);
         }
      } else {
         // b=90+(double)rank;
         b=0;
      }
      if(rank_down ==-1){j=ny+1; for (i=0; i<nx+2; i++){ U[i][j]=inity[i];};};
      if(rank_up   ==-1){j=0;    for (i=0; i<nx+2; i++){ U[i][j]=b;       };};
      if(rank_left ==-1){i=0;    for (j=0; j<ny+2; j++){ U[i][j]=b;       };};
      if(rank_right==-1){i=nx+1; for (j=0; j<ny+2; j++){ U[i][j]=b;       };};
   }
   //fflush(stdout);MPI_Barrier(MPI_COMM_WORLD);
   return 0;
}

/******************************************************************************/
/* euler  *****************************************/
int euler(){
   int i, j;
   int i_1, i_nx, j_1, j_ny;            // boundary indicies

   junk = laplacian();

   /*
   // if (i,j) not on a boundary
   if(rank_down ==-1){j_1 =2; } else {j_1 =1;   };
   if(rank_up   ==-1){j_ny=ny;} else {j_ny=ny+1;};
   if(rank_left ==-1){i_1 =2; } else {i_1 =1;   };
   if(rank_right==-1){i_nx=nx;} else {i_nx=nx+1;};

   //j_1 =1;
   //j_ny=ny+1;
   //i_1 =1;
   //i_nx=nx+1;
   */

   for (i=1; i<nx+1; i++){
      for (j=1; j<ny+1; j++){
         U[i][j] = U[i][j] + h*lu[i][j];
      }
   }

   return 0;
}

/******************************************************************************/
/* laplacian *****************************************/
int laplacian(){
   int i, j;
   junk = ghost_cell_update();
   if(sw_method==1){
      for (i=1; i<nx+1; i++){
         for (j=1; j<ny+1; j++){
            lu[i][j] = (U[i-1][j]+U[i+1][j]-4*U[i][j]+U[i][j-1]+U[i][j+1])/(pow(delta_x,2));
            //printf("i j lu... %3d %3d %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f %10.5f\n",i,j,lu[i][j],U[i-1][j],U[i+1][j],-4*U[i][j],U[i][j-1],U[i][j+1],pow(delta_x,2));
         }
      }
   }
   if(sw_method==2){
      for (i=1; i<nx+1; i++){
         for (j=1; j<ny+1; j++){
            // can speed up by raising to the power in own loop, then use those values here
            lu[i][j] = (pow(U[i-1][j],4)+pow(U[i+1][j],4)-4*pow(U[i][j],4)+pow(U[i][j-1],4)+pow(U[i][j+1],4))/(pow(delta_x,2));
         }
      }
   }
   return 0;
}

/******************************************************************************/
/* ghost_cell_update *****************************************/
int ghost_cell_update(){
   int i, j;
   int ierr, ierr1, ierr2, ierr3, tag=0;
   MPI_Status status_rreq[4], status_sreq[4];
   MPI_Request rreq[4], sreq[4];  // receive and send requests (down, up, left, right)
   double bufrd[nx], bufsd[nx], bufru[nx], bufsu[nx], bufrl[ny], bufsl[ny], bufrr[ny], bufsr[ny]; // receive and send buffers

   // DOWN **********************
   if(rank_down==-1){ // down boundary condition;
     //j=ny; for (i=1; i<=nx; i++){U[i][j]=U[i][j+1];}
      sreq[0] = MPI_SUCCESS; rreq[0] = MPI_SUCCESS;
   } else { // MPI
      j=ny; for (i=1; i<nx+1; i++){ bufsd[i-1]=U[i][j];}
      ierr2 = MPI_Isend(&bufsd, nx, MPI_DOUBLE, rank_down, tag, MPI_COMM_WORLD, &sreq[0]); if (ierr2 != MPI_SUCCESS) {printf("MPI_Isend ERROR status %d.\n", ierr2); ierr2=MPI_SUCCESS;}
      ierr1 = MPI_Irecv(&bufrd, nx, MPI_DOUBLE, rank_down, tag, MPI_COMM_WORLD, &rreq[0]); if (ierr1 != MPI_SUCCESS) {printf("MPI_Irecv ERROR status %d.\n", ierr1); ierr1=MPI_SUCCESS;}
   };

   // UP **********************
   if(rank_up==-1){ // up boundary condition;
     //j=1; for (i=1; i<=nx; i++){U[i][j]=U[i][j-1];}
      sreq[1] = MPI_SUCCESS; rreq[1] = MPI_SUCCESS;
   } else { // MPI
      j=1; for (i=1; i<nx+1; i++){ bufsu[i-1]=U[i][j];}
      ierr2 = MPI_Isend(&bufsu, nx, MPI_DOUBLE, rank_up, tag, MPI_COMM_WORLD, &sreq[1]); if (ierr2 != MPI_SUCCESS) {printf("MPI_Isend ERROR status %d.\n", ierr2); ierr2=MPI_SUCCESS;}
      ierr1 = MPI_Irecv(&bufru, nx, MPI_DOUBLE, rank_up, tag, MPI_COMM_WORLD, &rreq[1]); if (ierr1 != MPI_SUCCESS) {printf("MPI_Irecv ERROR status %d.\n", ierr1); ierr1=MPI_SUCCESS;}
   };

   // LEFT **********************
   if(rank_left==-1){ // left boundary condition;
     //i=1; for (j=1; j<=ny; j++){U[i][j]=U[i-1][j];}
      sreq[2] = MPI_SUCCESS; rreq[2] = MPI_SUCCESS;
   } else { // MPI
      i=1; for (j=1; j<ny+1; j++){ bufsl[j-1]=U[i][j];}
      ierr2 = MPI_Isend(&bufsl, ny, MPI_DOUBLE, rank_left, tag, MPI_COMM_WORLD, &sreq[2]); if (ierr2 != MPI_SUCCESS) {printf("MPI_Isend ERROR status %d.\n", ierr2); ierr2=MPI_SUCCESS;}
      ierr1 = MPI_Irecv(&bufrl, ny, MPI_DOUBLE, rank_left, tag, MPI_COMM_WORLD, &rreq[2]); if (ierr1 != MPI_SUCCESS) {printf("MPI_Irecv ERROR status %d.\n", ierr1); ierr1=MPI_SUCCESS;}
   };

   // RIGHT **********************
   if(rank_right==-1){ // right boundary condition;
     //i=nx; for (j=1; j<=ny; j++){U[i][j]=U[i+1][j];}
      sreq[3] = MPI_SUCCESS; rreq[3] = MPI_SUCCESS;
   } else { // MPI
      i=nx; for (j=1; j<ny+1; j++){ bufsr[j-1]=U[i][j];}
      ierr2 = MPI_Isend(&bufsr, ny, MPI_DOUBLE, rank_right, tag, MPI_COMM_WORLD, &sreq[3]); if (ierr2 != MPI_SUCCESS) {printf("MPI_Isend ERROR status %d.\n", ierr2); ierr2=MPI_SUCCESS;}
      ierr1 = MPI_Irecv(&bufrr, ny, MPI_DOUBLE, rank_right, tag, MPI_COMM_WORLD, &rreq[3]); if (ierr1 != MPI_SUCCESS) {printf("MPI_Irecv ERROR status %d.\n", ierr1); ierr1=MPI_SUCCESS;}
   };


   // update from receives
   ierr3 = MPI_Waitall(4, rreq, MPI_STATUSES_IGNORE); if (ierr3 != MPI_SUCCESS) {printf("MPI_Waitall ERROR status %d.\n", ierr3); ierr3=MPI_SUCCESS;}
   if(rank_down  != -1){j=ny+1;for (i=1; i<nx+1; i++){U[i][j]=bufrd[i-1];}};
   if(rank_up    != -1){j=0   ;for (i=1; i<nx+1; i++){U[i][j]=bufru[i-1];}};
   if(rank_left  != -1){i=0   ;for (j=1; j<ny+1; j++){U[i][j]=bufrl[j-1];}};
   if(rank_right != -1){i=nx+1;for (j=1; j<ny+1; j++){U[i][j]=bufrr[j-1];}};

   // finish waiting for sends
   ierr3 = MPI_Waitall(4, sreq, MPI_STATUSES_IGNORE); if (ierr3 != MPI_SUCCESS) {printf("MPI_Waitall ERROR status %d.\n", ierr3); ierr3=MPI_SUCCESS;}

   return 0;
}


int max_error(){
   int i, j;
   double diff, diff2, global_max;

   diff=0;
   for(j=1; j<ny+1; j++){
      for(i=1; i<nx+1; i++){
         diff2 = fabs(lu[i][j] - Ue[i][j]);
         if(diff<=diff2) diff = diff2;
         //printf("rank i j lu Ue diff2 %4d %4d %4d %14.6f %14.6f %14.6f\n",rank,i,j,lu[i][j],Ue[i][j],diff2);
      }
   }

   //printf("rank=%d local_max = %15.10f\n", rank, diff);

   MPI_Reduce(&diff, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);

   if(rank==0) printf("*** global_max = %15.10f\n", global_max);
   if(rank==0) printf("%dx%d %d %d %d %15.10f\n", nxg, nyg, nprocs, nprocsx, nprocsy, global_max);

   return 0;
}


int string2int(char *digit) {
   int result = 0;
   //--- Convert each digit char and add into result.
   while (*digit >= '0' && *digit <='9') {
      result = (result * 10) + (*digit - '0');
      digit++;
   }
   //--- Check that there were no non-digits at end.
   if (*digit != 0) {return -1;}
   return result;
}


/*
   output data from M471 heat equation problem

   input:   nx  number of interior points, x direction
            ny  number of interior points, y direction
            x[nx]  x coordinates of interior points
            y[ny]  y coordinates of interior points
            u[nx][ny]   solution at interiod points
            time

   Note: the inputs to this routine are vectors which only contain
   interior points - no ghost points.
   But this routine only accepts arrays of size u[nx][ny].

   For example, if you are following the convention I used in class,
   the u array will be something like:  u[nx+2][ny+2].
   where the interior points are:   u[2:nx+1][2:ny+1] and
   the ghost cells are u[1][:],u[nx+2][:], etc...


   so when calling this routine, in C you should first create new arrays:
   then copy the *interior* data only into these arrays, and then call this
   routine:

      double xoutput[nx],youtput[ny],uoutput[nx][ny]
      for (i=0; i<nx; ++i){
      for (j=0; j<ny; ++j){
           uoutput[i][j] = u[i+1][j+1];
      }
      }

      and similar for xoutput[] and youtput[].
      output(nx,ny,xoutput,youtput,uoutput,time)

   Call this routine once for each snapshot you wish to save to
   that it can be plotted in Matlab.

   The first time you call this routine, it will create the file
   and write (as IEEE double precision data) the input data:
   nx,ny,x,y,u,time

   subsequent calls will append to the same file
*/
int output(int nx,int ny,double x[],double y[],double u[nx][ny],double time) {

  double xout[nx],yout[ny],uout[nx][ny];
  double uout2[ny][nx];  // array where x varies fastest
  double tmp;
  int ierr,amode,rank,nproc,p,i,j;
  MPI_File fh;
  static int firstcall=1;
  char *fname="heat.out";

  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  MPI_Comm_size(MPI_COMM_WORLD,&nproc);

  if(rank==0) printf("Write: %s at t=%f\n",fname,time);

  if (firstcall) {
    amode = MPI_MODE_WRONLY | MPI_MODE_CREATE;
    MPI_File_delete(fname,MPI_INFO_NULL);
    firstcall=0;
  }else{
    amode = MPI_MODE_WRONLY | MPI_MODE_APPEND;
  }
  MPI_File_open(MPI_COMM_WORLD,fname,amode,MPI_INFO_NULL,&fh);

  if (rank==0) {
    double tmp=nproc;
    MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE);
    tmp=nx;
    MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE);
    tmp=ny;
    MPI_File_write(fh,&tmp,1,MPI_DOUBLE,MPI_STATUSES_IGNORE);
  }


  for (p=0; p<nproc; ++p) {
    if (rank==p) {
      if (p==0)  {
	memcpy(uout,u,8*nx*ny);
	memcpy(xout,x,8*nx);
	memcpy(yout,y,8*ny);
      }else{
	MPI_Send(x,nx,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
        MPI_Send(y,ny,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
        MPI_Send(u,nx*ny,MPI_DOUBLE,0,0,MPI_COMM_WORLD);
      }
    }
    if (rank==0) {
      if (p!=0) {
	MPI_Recv(xout,nx,MPI_DOUBLE,p,0,MPI_COMM_WORLD,MPI_STATUSES_IGNORE);
	MPI_Recv(yout,ny,MPI_DOUBLE,p,0,MPI_COMM_WORLD,MPI_STATUSES_IGNORE);
        MPI_Recv(uout,nx*ny,MPI_DOUBLE,p,0,MPI_COMM_WORLD,MPI_STATUSES_IGNORE);
      }
      // take the transpose of uout:
      for (i=0; i<nx; ++i) for (j=0; j<ny; ++j ) uout2[j][i]=uout[i][j];
      // output arrays:
      MPI_File_write(fh,xout,nx,MPI_DOUBLE,MPI_STATUSES_IGNORE);
      MPI_File_write(fh,yout,ny,MPI_DOUBLE,MPI_STATUSES_IGNORE);
      MPI_File_write(fh,uout2,nx*ny,MPI_DOUBLE,MPI_STATUSES_IGNORE);
    }
  }

  if (rank==0)
   MPI_File_write(fh,&time,1,MPI_DOUBLE,MPI_STATUSES_IGNORE);

  MPI_File_close(&fh);
  return 0;
}

/* eof */


/*

x=[0:.1:1];
xx=sin(pi*x)'*sin(pi*x);

yy=zeros(length(x));
for j=2:10;for i=2:10;
yy(i,j)=(xx(i-1,j)+xx(i+1,j)-4*xx(i,j)+xx(i,j-1)+xx(i,j+1))/.1^2;
end;end;

dif=yy- -2*pi^2*xx;

*/