/****************************************************************************
   * DESCRIPTION:  
   *   Serial HEAT2D Example - C Version
   *   This example is based on a simplified 
   *   two-dimensional heat equation domain decomposition.  The initial 
   *   temperature is computed to be high in the middle of the domain and 
   *   zero at the boundaries.  The boundaries are held at zero throughout 
   *   the simulation.  During the time-stepping, an array containing two 
   *   domains is used; these domains alternate between old data and new data.
   *
   *  The physical region, and the boundary conditions, are suggested
      by this diagram;
  
                     u = 0
               +------------------+
               |                  |
      u = 100  |                  | u = 100
               |                  |
               +------------------+
                     u = 100
  
  Interrior point :
    u[Central] = (1/4) * ( u[North] + u[South] + u[East] + u[West] )
  
   ****************************************************************************/
  #include <stdio.h>
  #include <stdlib.h>
  #include <math.h>
  #include <omp.h>
  
  #define N 500
  #define M 500  
  
  #define ITER_PRINT 100
  #define PRINT_DATA 1
  
  #define EPSILON 1e-1
  
  
  void update(int nx,int ny, float *u, float *unew, float * diff);
  void inidat(int nx, int ny, float *u, float *unew); 
  void prtdat(int nx, int ny, float *u,const char *fnam);
  
  
  
  
  int main(int argc, char *argv[])
  {
  
      float diff=1.0;
  
      float *u     = (float *)malloc(N * M * sizeof(float));
      float *unew  = (float *)malloc(N * M * sizeof(float));
   
      if(u==0 || unew ==0)
      {
          perror("Can't allocated data
  ");
          return -1;
      }
  
       printf ( "
  " );
       printf ( "HEATED_PLATE
  " );
       printf ( "  Parallel OpenMP version, using %d Threads
  ",omp_get_max_threads() );
       printf ( "  A program to solve for the steady state temperature distribution
  " );
       printf ( "  over a rectangular plate.
  " );
       printf ( "  Spatial grid of %d by %d points.
  
  ", M, N );
  
  
      /* Initialize grid and create input file */
      printf("Initializing grid
  ");
      
      inidat(N, M,u,unew);
  
      prtdat(N, M,u, "initial.dat");
      
     
      printf("Start computing
  ");
  
      int iter=0;
  
      /* 
       *   iterate until the  new solution unew differs from the old solution u
       *     by no more than EPSILON.
       *     */
       
      while(diff> EPSILON) {
  
          update(N, M, u, unew,&diff);
      
          if(iter%ITER_PRINT==0)
          printf("Iteration %d, diff = %f
   ", iter,diff);
  
          iter++;
      }
  
      prtdat(N, M, u, "final.dat");
       
      free(u);
      free(unew);
  }
  
  
  
  /****************************************************************************
   *  subroutine update
   ****************************************************************************/
  void update(int nx,int ny, float *u, float *unew, float * diff)
  {
      int ix, iy;
      *diff=0.0;
  
  #pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy)
      for (ix = 1; ix < nx-1; ix++) {
          for (iy = 1; iy < ny-1; iy++) {
              unew[ix*ny+iy] =  
                  (u[(ix+1)*ny+iy] +  u[(ix-1)*ny+iy] +
                   u[ix*ny+iy+1] +  u[ix*ny+iy-1] )/4.0;
            
          }
      }
  
  //compute reduction 
  
  
          float mydiff;
          
  #pragma omp parallel  shared(nx,ny,u,unew, diff) private (ix,iy,mydiff)
          {
  mydiff=0.0;
  #pragma omp for
      for (ix = 1; ix < nx-1; ix++) {
          for (iy = 1; iy < ny-1; iy++) {
              if (mydiff < fabs (unew[ix*ny+iy] - u[ix*ny+iy] ))
              {
                  mydiff = fabs ( unew[ix*ny+iy] - u[ix*ny+iy] );
              }
          }
      }
  
  
  # pragma omp critical
            {
                if (*diff < mydiff )
                {
                    *diff = mydiff;
                }
            }
  
  
  #pragma omp for
      for (ix = 1; ix < nx-1; ix++) {
          for (iy = 1; iy < ny-1; iy++) {
              u[ix*ny+iy] = unew[ix*ny+iy]; 
          }
      }   
  
  
  
  
          }
  }
  
  /*****************************************************************************
   *  Initialize Data
   *****************************************************************************/
  void inidat(int nx, int ny, float *u, float *unew) 
  {
      int ix, iy;
  
      /*
       *Set boundary data and interrior values
       * */
      for (ix = 0; ix < nx; ix++) 
          for (iy = 0; iy < ny; iy++) { 
  
              if(ix==0)
              {
                  u[ix*ny+iy]=0.0; 
              }
              else          
              if(iy==0 && ix!=0)
              {
                  u[ix*ny+iy]=100.0;
              }else
             
              if(ix==nx-1)
              {
                  u[ix*ny+iy]=100.0;
              }else
  
              if(iy==ny-1 && ix!=0)
              {   
                  u[ix*ny+iy]=100.0;
              }else
  
               u[ix*ny+iy]=( float ) ( 2 * nx + 2 * ny - 4 );
          }
  }
  
  /**************************************************************************
   * Print Data to files
   **************************************************************************/
  void prtdat(int nx, int ny, float *u,const char *fnam)
  {
   
      int ix, iy;
      FILE *fp;
  
      if(ITER_PRINT==0)return;
      
      fp = fopen(fnam, "w");
  
      for (ix = 0 ; ix < nx; ix++) {
          for (iy =0; iy < ny; iy++) {
  
              fprintf(fp, "%8.3f", u[ix*ny+iy]);
  
              if(iy!=ny-1)
              {
                  fprintf(fp, " ");
              }else
              {
                  fprintf(fp, "
  ");
              }
          }
      }
  
      fclose(fp);
  }