/****************************************************************************
 * DESCRIPTION:  
 *   OpenMP 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>

#define NN 50
#define MM 50  

#define ITER_PRINT 100
#define _MAX_ITER 400 
#define PRINT_DATA 1
#define _EPSILON 0.01


float update(int nx,int ny, float *u, float *unew);
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[])
{
    int N=NN,M=MM;
    float EPSILON=_EPSILON;
    int MAX_ITER= _MAX_ITER;

    if(argc !=4)
    {
        fprintf(stderr,"usage %s N EPSILON MAX_ITER\n ", argv[0]);
        fprintf(stderr,"\t\twhere N is GRID size, EPSILON is  Tolerance, MAX_ITER is max iteration\n"); 
        fprintf(stderr,"\t\texample N = 100, EPSILON = 0.1, MAX_ITER=1000\n"); 
        return -1;
    } 

    N = M = atoi(argv[1]);
    EPSILON = atof(argv[2]);
    MAX_ITER = atoi (argv[3]);

    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\n");
        return -1;
    }

    printf ( "\n" );
    printf ( "HEATED_PLATE\n" );
    printf ( "  OpenMP version\n" );
    printf ( "  A program to solve for the steady state temperature distribution\n" );
    printf ( "  over a rectangular plate.\n" );
    printf ( "\n" );
    printf ( "  Spatial grid of %d by %d points.\n", M, N );
    printf ( "  The iteration will end until tolerance <= %f\n\n",EPSILON);

    /* Initialize grid and create input file */
    printf("Initializing grid\n");

    inidat(N, M,u,unew);

    prtdat(N, M,u, "initial.dat");

    printf("Start computing\n\n");

    int iter=0;

    /* 
     *   iterate until the  new solution unew differs from the old solution u
     *     by no more than EPSILON.
     *     */

    while(diff> EPSILON && iter <MAX_ITER) {

        diff= update(N, M, u, unew);

        if(iter%ITER_PRINT==0)

            printf("Iteration %d, diff = %e\n ", iter,diff);

        iter++;
    }

    if(diff>EPSILON)
       printf("***Not converged MAX_ITERATION %d reached\n***", MAX_ITER);

    prtdat(N, M, u, "final.dat");

    free(u);
    free(unew);
}



/****************************************************************************
 *  subroutine update
 ****************************************************************************/
float update(int nx,int ny, float *u, float *unew)
{
    int ix, iy;
    float diff=0.0;

#pragma omp parallel for shared(nx,ny,u,unew) private (ix,iy) reduction(max:diff)
    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] )*0.25;
 diff = fmax( diff, fabs(unew[ix*ny+iy] - u[ix*ny+iy]));
        }
    }
   

/*
 * COPY OLD DATA
 */
#pragma omp parallel for
        for (ix = 1; ix < nx-1; ix++) {
            for (iy = 1; iy < ny-1; iy++) {
                u[ix*ny+iy] = unew[ix*ny+iy]; 
            }
        }   
    

    return diff;
}

/*****************************************************************************
 *  Initialize Data
 *****************************************************************************/

void inidat(int nx, int ny, float *u, float *unew) 
{
    int ix, iy;

    /*
     *Set boundary data and interrior values
     * */

#pragma omp parallel private (ix,iy)
    {
    // interior points
    #pragma omp for    
    for (ix = 1; ix < nx-1; ix++) 
        for (iy = 1; iy < ny-1; iy++) { 
            u[ix*ny+iy]=0.0; 
        }

   //boundary left
    #pragma omp for    
    for (ix = 1; ix < nx-1; ix++){ 
        u[ix*ny]=100.0; 

    }

    //boundary right
    #pragma omp for    
    for (ix = 1; ix < nx-1; ix++){ 
        u[ix*ny+ (ny-1)]=100.0; 

    }

    //boundary down
    #pragma omp for    
    for (iy = 0; iy < ny; iy++){ 
        u[(nx-1)*(ny)+iy]=100.0; 

    }

    //boundary top
    #pragma omp for 
    for (iy = 0; iy < ny; iy++){ 
        u[iy]=100.0; 

    }
   
    }

}





/**************************************************************************
 * 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, "\n");
            }
        }
    }

    fclose(fp);
}