#define NPES     4
#define NC       1000                   /* Number of Cols        */
#define NR       1000                   /* Number of Rows        */
#define NRL      NR/NPES                /* Number of Rows per PE */
#define DOWN     100                    /* Tag for messages down */
#define UP       101                    /* Tag for messages up   */
#define ROOT     0                      /* The root PE           */
#define MAX(x,y) ( ((x) > (y)) ? x : y )

#include <stdio.h>
#include <math.h>
#include <mpi.h>                        /* Required MPI library */

void initialize( float t[NRL+2][NC+2] );
void set_bcs   ( float t[NRL+2][NC+2], int mype, int npes );

int main( int argc, char **argv ){
  
  int        npes;                      /* Number of PEs */
  int        mype;                      /* My PE number  */
  int        stat;                      /* Error Status  */
  int        niter;                     /* iter counter  */
  MPI_Status status;

  float      t[NRL+2][NC+2], told[NRL+2][NC+2];
  float      dt;                        /* Delta t       */
  float      dtg;                       /* Delta t global*/

  int        i, j, iter;

   /*-----------------*/
  /*  Initialize MPI */
 /*-----------------*/

  MPI_Init(&argc, &argv);

   /*----------------------------------------*/
  /*  Determine size of global communicator */
 /*----------------------------------------*/

  MPI_Comm_size(MPI_COMM_WORLD, &npes );

   /*-----------------------------------------------*/
  /*  Determine my rank in the global communicator */
 /*-----------------------------------------------*/

  MPI_Comm_rank(MPI_COMM_WORLD, &mype );

  initialize(t);                  /* Give initial guesss of 0. */

  set_bcs(t, mype, npes);         /* Set the Boundary values   */

  for( i=0; i<=NRL+1; i++ )       /* Copy the values into told */
    for( j=0; j<=NC+1; j++ )
      told[i][j] = t[i][j];

   /*-------------------------------------------------*/
  /* Do Computation on Sub-grid for Niter iterations */
 /*-------------------------------------------------*/

  niter=1000;

  for( iter=1; iter<=niter; iter++ ) {

    for( i=1; i<=NRL; i++ )
      for( j=1; j<=NC; j++ )
        t[i][j] = 0.25 * ( told[i+1][j] + told[i-1][j] +
                           told[i][j+1] + told[i][j-1] );
    dt = 0.;

    for( i=1; i<=NRL; i++ )       /* Copy for next iteration  */
      for( j=1; j<=NC; j++ ){
	dt         = MAX( abs(t[i][j]-told[i][j]), dt);
        told[i][j] = t[i][j];
      }

    if( mype < npes-1 )

   /*------------------------------------------------------------------*/
  /* Send my data down to the processor below me; Only npes-1 do this */
 /*------------------------------------------------------------------*/

      MPI_Send(&t[NRL][1], NC, MPI_FLOAT, mype+1, DOWN, MPI_COMM_WORLD);

    if( mype != 0 )

   /*--------------------------------------------------------------------*/
  /* Sending my data up to the processor above me ; Only npes-1 do this */
 /*--------------------------------------------------------------------*/

      MPI_Send(&t[1][1], NC, MPI_FLOAT, mype-1, UP, MPI_COMM_WORLD);

    if( mype != 0 )

   /*---------------------------------------------------*/
  /*  Receive new data from UP processor of any source */
 /*---------------------------------------------------*/

      MPI_Recv(&t[0][1], NC, MPI_FLOAT, MPI_ANY_SOURCE, DOWN,
               MPI_COMM_WORLD, &status);

    if( mype != npes-1 )

   /*------------------------------------------------------*/
  /*  Receive new data from DOWN processsor of any source */
 /*------------------------------------------------------*/

      MPI_Recv(&t[NRL+1][1], NC, MPI_FLOAT, MPI_ANY_SOURCE, UP,
               MPI_COMM_WORLD, &status);

   /*---------------------------------------------------------------------*/
  /*  Detemermine max value from all processor calculations using reduce */
 /*---------------------------------------------------------------------*/

      MPI_Reduce(&dt, &dtg, 1, MPI_FLOAT, MPI_MAX, ROOT, MPI_COMM_WORLD);

   /*------------------------*/
  /* Print some test values */
 /*------------------------*/

    if( (iter%100) == 0 ) {
      if( mype == 0 )
        printf("Iter = %4d: PE = %d: t[10][10] = %20.8f\n", 
                iter, mype, t[10][10]); 
    }  

   /*---------------------------------------------------------------*/
  /* All processors in the global communicator wait at the barrier */
 /*---------------------------------------------------------------*/

    MPI_Barrier( MPI_COMM_WORLD );

  }  /* End of iteration */

  MPI_Finalize();       /* Finalize MPI */

}    /* End of Program */

   /*-----------------------------------------------------*/
  /* Initialize all the values to 0. as a starting value */
 /*-----------------------------------------------------*/

void initialize( float t[NRL+2][NC+2] ){

  int        i, j, iter;
	  
  for( i=0; i<=NRL+1; i++ )        /* Initialize */
    for ( j=0; j<=NC+1; j++ )
      t[i][j] = 0.0;
}

    /*----------------------------------------------------------------*/
   /* Set the values at the boundary.  Values at the boundary do not */
  /* Change through out the execution of the program                */
 /*----------------------------------------------------------------*/

void set_bcs( float t[NRL+2][NC+2], int mype, int npes ){

  int i, j;

  for( i=0; i<=NRL+1; i++ ) {      /* Set Left and Right boundary */
    t[i][0]       = 100.0;
    t[i][NC+1] = 100.0;
  }

  if( mype == 0 )                  /* Top boundary */
    for( j=0; j<=NC+1; j++ )
      t[0][j]        = 100.0;

  if( mype == npes-1 )             /* Bottom boundary */
    for( j=0; j<=NC+1; j++ )
      t[NRL+1][j] = 100.0;

}