#include <stdio.h>
#include "mpi.h"
#define N 100
#define Tol 0.00001

float **matrix(int m, int n)
{  int i;
   float **ptr;
   ptr = (float **)calloc(m,sizeof(float *));
   for(i=0; i < m; i++)
         ptr[i] = (float *) calloc(n, sizeof(float));
   return (ptr);
}

float iteration( float **old, float **new, int start, int finish)
{   float diff, maxerr =0;
     int i,j;
     for (i= start; i< finish; i++)
          for (j=1; j<N; j++){ new[i][j] = 0.25*(old[i+1][j] + old[i-1][j] + 
                                                   old[i][j+1] +old[i][j-1]);
     diff = new[i][j] - old[i][j];
                    if(diff < 0) diff = -diff;
                    if(maxerr < diff) maxerr= diff;      }
     return (maxerr);
}

main (int argc, char **argv)
{     
      float **new, **old, **tmp, maxerr, err, maxerrG;
      int nprocs, nid, remainder, size, i, j;
      char str[20];
      FILE *fp;
      MPI_Status status;
      MPI_Request req_send10, req_send20, req_recv10, req_recv20;
      MPI_Init(&argc, &argv);
      MPI_Comm_rank(MPI_COMM_WORLD, &nid);
      MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
      remainder = (N-1)% nprocs;
      size = (N-1-remainder)/nprocs;
      if(nid < remainder)
             size=size+2;
      else
             size = size+1; 
      new = matrix (size+1, N+1);
      old  = matrix (size+1, N+1);
      for (i=0; i < size +1; i++)
          new[i][0] = new [i][N] = old[i][0] = old [i][N] = 1.0;
      if(nid == 0)
      for(j=1; j< N; j++)
          new[0][j] = old[0][j]  =1.0;
      if(nid == nprocs-1)
          for (j=1; j<N ; j++)
                new[size][j] = old[size][j] =1.0;
      maxerr = iteration (old, new, 1, size);
      MPI_Allreduce(&maxerr, &maxerrG, 1, MPI_FLOAT, MPI_MAX, MPI_COMM_WORLD);
      while (maxerrG > Tol){
          tmp = new;
          new = old;
          old = tmp;
          req_send10 = req_recv20 = MPI_REQUEST_NULL;
          if(nid < nprocs-1){
               MPI_Isend(&old[size-1][1], N-1,MPI_FLOAT, nid+1,10, 
                                  MPI_COMM_WORLD, &req_send10);
               MPI_Irecv(&old[size][1], N-1,MPI_FLOAT, nid+1, 20, 
                                  MPI_COMM_WORLD, &req_recv20);
                   }
          req_send20 = req_recv10 = MPI_REQUEST_NULL;
          if (nid > 0) {
               MPI_Isend(&old[1][1], N-1, MPI_FLOAT, nid-1, 20, 
                                 MPI_COMM_WORLD, &req_send20);
               MPI_Irecv(&old[0][1], N-1, MPI_FLOAT, nid-1, 10, 
                                 MPI_COMM_WORLD, &req_recv10);
                 }
         maxerr = iteration (old, new, 2, size-1);
         if(nid < nprocs-1)
                  MPI_Wait(&req_recv20,&status);
         err = iteration (old,new,size-1,size);
         if(err > maxerr) maxerr = err;
         if(nid >0)
                  MPI_Wait(&req_recv10,&status);
         err = iteration (old,new,1,2);
         if(err > maxerr) maxerr = err;
         MPI_Allreduce(&maxerr,&maxerrG, 1,MPI_FLOAT,MPI_MAX, MPI_COMM_WORLD);
     }
      sprintf (str, "Solution%d.Txt",nid);
      fp = fopen(str,"wt");
 
if(nid ==0)
     for(j=0; j <N+1; j++)
           fprintf(fp,"%6.4f\n",new[0][j]);
for(i=1; i < size; i++)
     for(j=0; j<N+1; j++)
            fprintf(fp, "%6.4f\n",new[i][j]);
if(nid == nprocs-1)
     for(j=0; j<N+1; j++)
            fprintf(fp, "%6.4f\n", new[size][j]);
fclose(fp);
MPI_Finalize();
}