#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 noprocs, 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_FLOAT, &nid);
    MPI_Comm_size(MPI_FLOAT, &noprocs);
    remainder = (N - 1) % noprocs;
    size = (N - 1 - remainder)/noprocs;
    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;
    if(nid == 0)
       for(j = 1; j < N; j++)
          new[0][j] = old[0][j] = 1;
    if(nid == noprocs - 1)
       for(j = 1; j < N; j++)
          new[size][j] = old[size][j] = 1;
    maxerr = iteration(old,new,1,size);
    MPI_Allreduce(&maxerr,&maxerrG,1,MPI_MAX,
                  MPI_COMM_WORLD);
    while(maxerrG > Tol){
       tmp = new;
       new = old;
       old = tmp;
       req_send10 = req_recv20 = MPI_REQUEST_NULL;
       if(nid < noprocs-1){
          MPI_Isnd(&old[size-1][1],N-1,MPI_FLOAT,nid+1,10,
                    MPI_COMM_WORLD,&req_send10);
          MPI_Ircv(&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_Isnd(&old[1][1],N-1,MPI_FLOAT,nid-1,20,
                    MPI_COMM_WORLD,&req_send20);
          MPI_Ircv(&old[0][1],N-1,MPI_FLOAT,nid-1,10,
                    MPI_COMM_WORLD,&req_recv10);
       }
       maxerr = iteration(old,new,2,size-1);
       if(nid < noprocs-1)
          MPI_Wait(MPI_COMM_WORLD,&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);
    }
    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 == noprocs - 1)
       for(j = 0; j < N + 1; j++)
          fprintf(fp,"%6.4f\n",new[size][j]);
    fclose(fp);
}