#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>

/* global objects: */

int* iters;
double XINC   = 0.0, YINC   = 0.0;
double Sreal  = 0.0, Simag  = 0.0;
int XYdots  =  0;
int NITER  = 1000;
int Procs = 1;



int MINVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v > a[e] ) v = a[e];
   
   }
   
   return(v);
}

int MAXVAL(int s, int*a)
{
   int v, e;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v < a[e] ) v = a[e];
   
   }
   
   return(v);
}
   
   
void Terminate(int p)
{
#include "mpi.h"
   int ierr;
   
   fprintf(stdout,"Process %d  terminates\n",p);
   ierr = MPI_Finalize();
   
   return;
}
   

void calcola( int X0,  int X1,  int Y0,  int Y1)
{

   int pos;
   int       IX, IY, IZ;
   double     CA, CB;
   double     ZA, ZB;
   double     ZSIZE ;
   double     ZNEWA, ZNEWB;

   for ( IX = X0; IX < X1; IX++ ) {
     /* New column  */
      CA = XINC * IX + Sreal ;
      for ( IY = Y0; IY < Y1; IY++ ) {
         ZA = 0.0  ;
         ZB = 0.0  ;
         CB = YINC * IY + Simag  ;
         ZSIZE = (CA*CA + CB*CB )  ;
         /* Iteration loop for each ( Z**2 + C )   */
         for ( IZ = 0; IZ <= NITER ; IZ++ ) {
            if ( ZSIZE >= 4 ) break; 
            ZNEWA = ( ZA*ZA ) - ( ZB*ZB ) + CA  ;
            ZNEWB = ( 2*ZA*ZB ) + CB  ;
            ZSIZE = ( ZNEWA*ZNEWA + ZNEWB*ZNEWB )  ;
            ZA   = ZNEWA  ;
            ZB   = ZNEWB  ;
         }  
         /* Results  */
         iters[IY*XYdots+IX] = IZ;
      } 
   }  
}

void IntData2pgm(int s1, int s2, int* idata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PGM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs;
   float  rmin, rmax;
   char  fname[80];
   /*
      Write on unit 700 with  PGM format
   */
   strcpy(fname,name);
   strcat(fname,".pgm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P2\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   fprintf(stdout,"rmin, rmax = %f %f\n",rmin, rmax);
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         vp = (int) ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vs++;
         fprintf(ouni,"%4d ", vp);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

void IntData2ppm(int s1, int s2, int* idata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PPM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs, c, r, g, b;
   float  fvp, rmin, rmax;
   char  fname[80];

   strcpy(fname,name);
   strcat(fname,".ppm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P3\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,idata); rmax = MAXVAL(s1*s2,idata);
   fprintf(stdout,"rmin, rmax = %f %f\n",rmin, rmax);
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         fvp = ( (idata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vp = (int)fvp % 256;
         c = vp / 42;   /* colour */
         switch (c) {
            case 0:
               vp = (int) (fvp - 1.0) * 255.0 / 42.0;
               r = vp; g = 0; b = 0;
               break;
            case 1:
               vp = (int) (fvp - 42.0) * 255.0 / 42.0;
               r = vp; g = vp; b = 0;
               break;
            case 2:
               vp = (int) (fvp - 84.0) * 255.0 / 42.0;
               r = 0; g = vp; b = 0;
               break;
            case 3:
               vp = (int) (fvp - 126.0) * 255.0 / 42.0;
               r = 0; g = vp; b = vp;
               break;
            case 4:
               vp = (int) (fvp - 168.0) * 255.0 / 42.0;
               r = 0; g = 0; b = vp;
               break;
            case 5:
               vp = (int) (fvp - 210.0) * 255.0 / 42.0;
               r = vp; g = 0; b = vp;
               break;
            case 6:
               r = 0; g = 0; b = fvp/10.0;
               break;
            default:
               r = 0; g = 0; b = 0;
         }
         vs++;
         fprintf(ouni,"%4d %4d %4d ", r,g,b);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}


/*
   vs = 0
   do i = 1, s1
      do j = 1, s2
         rvp = ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         ! Palette: r,g,b values
         vp = MOD(INT(rvp),256)
         SELECT CASE (vp)
            CASE(1:39)
               vp = INT ( (rvp - 1.0) * 255.0 / 40.0 )
               r = vp; g = 0; b = 0
            CASE(40:79)
               vp = INT ( (rvp - 40.0) * 255.0 / 40.0 )
               r = vp; g = vp; b = 0
            CASE(80:119)
               vp = INT ( (rvp - 80.0) * 255.0 / 40.0 )
               r = 0; g = vp; b = 0
            CASE(120:159)
               vp = INT ( (rvp - 120.0) * 255.0 / 40.0 )
               r = 0; g = vp; b = vp
            CASE(160:199)
               vp = INT ( (rvp - 160.0) * 255.0 / 40.0 )
               r = 0; g = 0; b = vp
            CASE(200:239)
               vp = INT ( (rvp - 200.0) * 255.0 / 40.0 )
               r = vp; g = 0; b = vp
            CASE(240:255)
               r = 0; g = 0; b = (rvp/10.0)
            CASE DEFAULT    ! 0 value
               r = 0; g = 0; b = 0
         END SELECT
         vs = vs + 1
         write(ouni,"(3(i4,1x))",advance="no") r, g, b
         if (vs >= 10 ) {
            write(ouni,"(a)") " "
            vs = 0
         }
      }
      write(ouni,"(a)") " "
      vs = 0
   }
   close(ouni)
   return
end subroutine IntData2ppm


END MODULE global

*/
/*

   MAND - Compute the points outside the Mandelbrot set
 
    Computes the iteration number for each pixel in a given squared area
 
    Syntax: 
    mandel Sreal Simag XYdots Range procs NITER 
               
  INPUT:            
     REAL(8) Sreal = Real coord. of lower left vertex 
     REAL(8) Simag = Imaginary coord. of lower left vertex 
     INTEGER XYdots = Horizontal and vertical pixels 
     INTEGER Range = Square area side length 
     INTEGER procs = Independent computing units 
     INTEGER NITER = Maximum iterations
               
 ************************************************************************** */

int main( int argc, char *argv[])    /* Mandel */

{
#include "mpi.h"

   /*  Mandel variables */
   double Range = 0.0;  
   double tempo2, tempo1, tempro;
   char  fname[40];

   int dX, dY, t, X0, X1, Y0, Y1;
   int st;
   FILE* inp_unit;
   
   /*  MPI variables */
   char server[MPI_MAX_PROCESSOR_NAME];
   int ierr, my_rank, numprocs, ls, dest, tag;
   MPI_Status status;
   int LetsGo, rc;

   ierr = MPI_Init(&argc,&argv);
   ierr = MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
   ierr = MPI_Comm_size(MPI_COMM_WORLD,&numprocs);
   ierr = MPI_Get_processor_name(server,&ls);

   LetsGo = 1;
   if ( my_rank == 0 ) {
   	   /* Master process looks for input data */
      tempo1 = MPI_Wtime();

      strcpy(fname,"mandel.inp");
      inp_unit = fopen(fname,"r");
      if ( ! inp_unit ) {
         fprintf(stderr,"!!!! Error read access to file %s\n",fname);
      }

      ierr = fscanf(inp_unit,"%lf, %lf, %d, %lf, %d, %d\n",&Sreal, &Simag, &XYdots, &Range, &Procs, &NITER);
      fprintf(stdout,"%lf, %lf, %d, %lf, %d, %d\n",Sreal, Simag, XYdots, Range, Procs, NITER);
      fprintf(stdout," ierr = %d\n",ierr);
/*
      if ( ierr != 0 ) {
   	     fprintf(stderr, "Error reading input data. Examples are:\n");
   	     fprintf(stderr, "-2.000, -2.000, 1000, 4.000, 1, 1000\n" );
   	     fprintf(stderr, "-1.500, -0.500, 1000, 0.500, 1, 1000\n" );
   	     fprintf(stderr, "-1.250, -0.200, 1000, 0.100, 1, 1000\n" );
   	     fprintf(stderr, "-1.250, -0.180, 1000, 0.030, 1, 1000\n" );
   	     fprintf(stderr, "-1.235, -0.171, 1000, 0.006, 1, 1000\n" );
   	     fprintf(stderr, "-1.233, -0.169, 1000, 0.001, 1, 1000\n" );
   	     fprintf(stderr, "-1.233, -0.169, 1000, 0.0001, 1, 1000\n" );
   	     LetsGo = 0;
      }
*/
      if ( LetsGo ) {
         /*  Master allocates entire matrix */
         iters = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating ITERS(%d,%d)\n",XYdots,XYdots);
      	     LetsGo = 0;
         }
      }
   }
   /* Every process synchronizes */
   ierr = MPI_Barrier(MPI_COMM_WORLD);
   /* Send start flag */
   ierr = MPI_Bcast(&LetsGo,1,MPI_INT,0,MPI_COMM_WORLD);
   if ( ierr != 0 ) {
      fprintf(stderr, " Error MPI_Bcast process %d \n",my_rank);
   }
   if ( ! LetsGo ) Terminate(my_rank);
      /* Send data to processes: */
      /*	  Sreal, Simag, XYdots, Range, NITER */
   ierr = MPI_Bcast(&Sreal,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&Simag,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&XYdots,1,MPI_INT,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&Range,1,MPI_DOUBLE,0,MPI_COMM_WORLD);
   ierr = MPI_Bcast(&NITER,1,MPI_INT,0,MPI_COMM_WORLD);
   fprintf(stderr, "Running MPI Mandel on server %s . Input data are: %lf %lf %d %lf %d %d\n", 
              server, Sreal, Simag, XYdots, Range, Procs, NITER);

   /* Every process executes */
   Procs = numprocs;   /* procs data required for compatibility reasons */

   /* ---> Compute complex spacing between pixels <---  */
   XINC = Range / (float)XYdots;
   YINC = Range / (float)XYdots;
  
/*  Let's decompose on 2nd dimension  */
/*   dX = ( XYdots + (Procs - 1) ) / Procs */
     dY = ( XYdots + (Procs - 1) ) / Procs;
       
       X0 = 0; X1 = XYdots;

       Y0 = dY * my_rank;
       Y1 = dY * ( my_rank + 1);
       if ( Y1 > XYdots ) Y1 = XYdots;
        
       if ( my_rank != 0 ) {   	   
         /*  Processes allocate local matrix */
         iters = (int*) malloc(4*XYdots*XYdots);
         if ( iters == NULL ) {
   	        fprintf(stderr, "Error allocating ITERS(%d,%d)\n",XYdots,XYdots);
                ! MPI_Abort should terminate all processes
                ierr = MPI_Abort(MPI_COMM_WORLD,rc);
                return(1);
         }
       }
       calcola(X0, X1, Y0, Y1);
       if ( my_rank != 0 ) {
          /*   Send data to master */
          dest = 0; tag = 100;
          ierr = MPI_Send(&iters[Y0*XYdots], (XYdots*(Y1-Y0)), MPI_INT, dest, tag, MPI_COMM_WORLD);
       }
	
   if ( my_rank == 0 ) {
      /* Master process gathers results */
      dX = ( XYdots + (Procs - 1) ) / Procs;
      dY = ( XYdots + (Procs - 1) ) / Procs;

      for ( t = 1; t < Procs; t++ ) {
        X0 = 0; X1 = XYdots;

        Y0 = dY * t;
        Y1 = dY * ( t + 1);
        if ( Y1 > XYdots ) Y1 = XYdots;
        tag = 100;
        ierr = MPI_Recv(&iters[Y0*XYdots], (XYdots*(Y1-Y0)), MPI_INT, t, tag, MPI_COMM_WORLD, &status);
      }
   
      tempo2 = MPI_Wtime(); 
      tempro =  (tempo2 - tempo1);

      IntData2pgm(XYdots,XYdots,iters,"mandel");
      IntData2ppm(XYdots,XYdots,iters,"mandel");
   
      fprintf(stdout, "Elapsed time is %lf \n",tempro);
   	   	   
   }
   MPI_Barrier(MPI_COMM_WORLD);
   Terminate(my_rank);
   return(0);
}