MODULE global
   IMPLICIT NONE
   INTEGER, DIMENSION(:,:), ALLOCATABLE :: iters

!  global objects: 
   REAL(8) :: XINC   = 0.0 !  Distance between pixels     
   REAL(8) :: YINC   = 0.0 !  Distance between pixels     
   REAL(8) :: Sreal  = 0.0 !  Beginning coordinate. Real.     
   REAL(8) :: Simag  = 0.0 !  Beginning coordinate. Imaginary.  
   INTEGER :: XYdots  =  0  ! Number of display pixels in X, Y directions     
   INTEGER :: NITER  = 1000!  Maximum number of iterations
   INTEGER :: Procs = 1
   
CONTAINS

subroutine IntData2pgm(s1,s2,idata,name)
!  Simple subroutine to dump integer data in a PGM format
   implicit none
   INTEGER, INTENT(IN) :: s1, s2
   INTEGER, DIMENSION(s1,s2), intent(in) :: idata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni
   integer :: vp, vs
   real :: rmin, rmax
   character(80) :: fname
!
!  Write on unit 700 with  PGM format
   ouni = 700
   fname = TRIM(name)//".pgm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      print*,"!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P2"
!  Dimensions
   write(ouni,*) s1, s2
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(idata); rmax = MAXVAL(idata)
   vs = 0
   do i = 1, s1
      do j = 1, s2
         vp = INT ( (idata(i,j) - rmin) * 255.0 / (rmax - rmin) )
         vs = vs + 1
         write(ouni,"(i4)",advance="no") vp
         if (vs >= 10 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
   end do
   close(ouni)
   return
end subroutine IntData2pgm


subroutine IntData2ppm(s1,s2,idata,name)
!  Simple subroutine to dump integer data in a PPM format
   implicit none
   INTEGER, INTENT(IN) :: s1, s2
   INTEGER, DIMENSION(s1,s2), intent(in) :: idata
   character(*), intent(in) :: name
!
   integer :: i, j, ouni
   integer :: vp, vs, r, g, b
   real :: rmin, rmax, rvp
   character(80) :: fname
!
!  Write on unit 700 with  PPM format
   ouni = 700
   fname = TRIM(name)//".ppm"
   open(ouni,file=trim(fname),status="replace",iostat=vs)
   if ( vs /= 0 ) then
      print*,"!!!! Error write access to file ",trim(fname)
   end if
!  Magic code
   write(ouni,"(a)") "P3"
!  Dimensions
   write(ouni,*) s1, s2
!  Maximum value
   write(ouni,"(a)") "255"
!  Values from 0 to 255
   rmin = MINVAL(idata); rmax = MAXVAL(idata)
   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 ) then
            write(ouni,"(a)") " "
            vs = 0
         end if
      end do
      write(ouni,"(a)") " "
      vs = 0
   end do
   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
!               
! ************************************************************************** 

PROGRAM Mandel 
   USE global
   IMPLICIT NONE
   INCLUDE "mpif.h"

!  Mandel variables
   REAL(8) :: Range=0.0, tempo2, tempo1, tempro
   CHARACTER(40) :: fname

   INTEGER :: dX, dY, t, X0, X1, Y0, Y1
   INTEGER :: st
   
!  MPI variables
   CHARACTER(MPI_MAX_PROCESSOR_NAME) :: server
   CHARACTER(32) :: params
   INTEGER :: ierr, my_rank, numprocs, ls, dest, tag, pos
   INTEGER status(MPI_STATUS_SIZE)
   LOGICAL :: LetsGo
!
   call MPI_INIT( ierr )
   call MPI_COMM_RANK( MPI_COMM_WORLD, my_rank, ierr )
   call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, ierr )
   call MPI_Get_processor_name(server,ls,ierr)
!
   LetsGo = .TRUE.
   IF ( my_rank == 0 ) THEN
   	   ! Master process looks for input data
      tempo1 = MPI_Wtime(ierr)

      fname = "mandel.inp"
      OPEN(11,FILE=fname,STATUS="OLD", IOSTAT=st)
      IF ( st /= 0 ) THEN
   	     WRITE(*,*) "Error access file ",TRIM(fname)
   	     LetsGo = .FALSE.
      END IF

      READ(11,*,IOSTAT=st) Sreal, Simag, XYdots, Range, procs, NITER
      IF ( st /= 0 ) THEN
   	     WRITE(*,*) "Error reading input data. Examples are:"
   	     WRITE(*,*) "-2.000, -2.000, 1000, 4.000, 1, 1000" 
   	     WRITE(*,*) "-1.500, -0.500, 1000, 0.500, 1, 1000" 
   	     WRITE(*,*) "-1.250, -0.200, 1000, 0.100, 1, 1000" 
   	     WRITE(*,*) "-1.250, -0.180, 1000, 0.030, 1, 1000" 
   	     WRITE(*,*) "-1.235, -0.171, 1000, 0.006, 1, 1000" 
   	     WRITE(*,*) "-1.233, -0.169, 1000, 0.001, 1, 1000" 
   	     WRITE(*,*) "-1.233, -0.169, 1000, 0.0001, 1, 1000" 
   	     LetsGo = .FALSE.
      END IF
      IF ( LetsGo ) THEN
         !  Master allocates entire matrix
         ALLOCATE(iters(0:XYdots-1,0:XYdots-1),STAT=st)
         IF ( st /= 0 ) THEN
   	        WRITE(*,*) "Error allocating ITERS(",XYdots,",",XYdots,")"
      	     LetsGo = .FALSE.
         END IF
      END IF
   END IF
   ! Every process synchronizes
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   ! Send start flag
   call MPI_BCAST(LetsGo,1,MPI_LOGICAL,0,MPI_COMM_WORLD,ierr)
   IF ( .NOT.LetsGo ) CALL Terminate(my_rank)
      ! Send data to processes:
      !	  Sreal, Simag, XYdots, Range, NITER
   IF ( my_rank == 0 ) THEN
      pos=0
      CALL MPI_Pack(Sreal,1,MPI_DOUBLE_PRECISION,params,32,pos,MPI_COMM_WORLD,ierr)
      CALL MPI_Pack(Simag,1,MPI_DOUBLE_PRECISION,params,32,pos,MPI_COMM_WORLD,ierr)
      CALL MPI_Pack(XYdots,1,MPI_INTEGER,params,32,pos,MPI_COMM_WORLD,ierr)
      CALL MPI_Pack(Range,1,MPI_DOUBLE_PRECISION,params,32,pos,MPI_COMM_WORLD,ierr)
      CALL MPI_Pack(NITER,1,MPI_INTEGER,params,32,pos,MPI_COMM_WORLD,ierr)
   END IF
   call MPI_BCAST(params,32,MPI_PACKED,0,MPI_COMM_WORLD,ierr)
   IF ( my_rank /= 0 ) THEN
      pos=0
      CALL MPI_Unpack(params,32,pos,Sreal,1,MPI_DOUBLE,MPI_COMM_WORLD,ierr)
      CALL MPI_Unpack(params,32,pos,Simag,1,MPI_DOUBLE,MPI_COMM_WORLD,ierr)
      CALL MPI_Unpack(params,32,pos,XYdots,1,MPI_INT,MPI_COMM_WORLD,ierr)
      CALL MPI_Unpack(params,32,pos,Range,1,MPI_DOUBLE,MPI_COMM_WORLD,ierr)
      CALL MPI_Unpack(params,32,pos,NITER,1,MPI_INT,MPI_COMM_WORLD,ierr)
   
   ELSE    
      WRITE(*,"(A,I4,2A)") "Running ",numprocs," 'MPI Mandel' processors on server ",TRIM(server)
      WRITE(*,"(A,2F12.8,I6,F12.4,I4,I8)") " Input data are: ", &
   &           Sreal, Simag, XYdots, Range, procs, NITER
   END IF
   ! Every process executes
   procs = numprocs

   ! ---> Compute complex spacing between pixels <--- 
   XINC = Range / real(XYdots)  
   YINC = Range / real(XYdots)  
  
!  Let's decompose on 2nd dimension 

   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 ) THEN   	   
         ALLOCATE(iters(0:XYdots-1,0:XYdots-1),STAT=st)
       END IF
       CALL calcola(X0, X1, Y0, Y1)
       IF ( my_rank /= 0 ) THEN
       !   Send data to master
          dest = 0; tag = 100
          call MPI_Send(iters(0,Y0), (XYdots*(Y1-Y0)), MPI_INTEGER, dest, &
                       tag, MPI_COMM_WORLD, ierr)
       END IF
	
   IF ( my_rank == 0 ) THEN
   	   ! Master process gathers results
   	  dX = ( XYdots + (Procs - 1) ) / Procs
      dY = ( XYdots + (Procs - 1) ) / Procs

     DO t = 1, Procs-1
       X0 = 0; X1 = XYdots

       Y0 = dY * t
       Y1 = dY * ( t + 1)
       if ( Y1 > XYdots ) Y1 = XYdots
       tag = 100
       call MPI_Recv(iters(0,Y0), (XYdots*(Y1-Y0)), MPI_INTEGER, t, &
                       tag, MPI_COMM_WORLD, status, ierr)       
      END DO
   
      tempo2 = MPI_Wtime(ierr)
   
      tempro =  (tempo2 - tempo1)

      CALL IntData2pgm(XYdots,XYdots,iters,"mandel")
      CALL IntData2ppm(XYdots,XYdots,iters,"mandel")
   
      WRITE(*,*) "Elapsed time is ",tempro
   	   	   
   END IF
   CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
   CALL Terminate(my_rank)
   STOP
END PROGRAM MANDEL  

SUBROUTINE Terminate(p)
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: P
   INCLUDE "mpif.h"
   INTEGER :: ierr
   
   WRITE(*,*) "Process ",p," terminates"
   CALL MPI_Finalize(ierr)
   
   STOP
END SUBROUTINE Terminate
   

SUBROUTINE calcola( X0,  X1,  Y0,  Y1)
   USE global
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: x0, x1, y0, y1 

   INTEGER :: ix, iy, iz
   REAL(8) :: ca, cb, za, zb
   REAL(8) :: zsize, znewa, znewb

   DO ix = x0, x1-1 
     ! New column 
      ca = XINC * ix + Sreal 
      DO iy = y0, y1-1 
         za = 0.0  
         zb = 0.0  
         cb = YINC * iy + Simag  
         zsize = (ca*ca + cb*cb )  
         ! Iteration loop for each ( Z**2 + C )  
         DO iz = 0, NITER
            IF ( zsize >= 4 ) EXIT 
            znewa = ( za*za ) - ( zb*zb ) + ca  
            znewb = ( 2*za*zb ) + cb  
            zsize = ( znewa*znewa + znewb*znewb )  
            za   = znewa  
            zb   = znewb  
         END DO  
         ! Results 
         iters(ix,iy) = iz
      END DO 
   END DO  
END SUBROUTINE calcola