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)
   print*,"rmin, rmax = ", rmin, rmax
   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)
   print*,"rmin, rmax = ", rmin, rmax
   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  !  Zoom factor       
   REAL ::   tempo2, tempo1, cpu2, cpu1, tempro, cpupro
   CHARACTER(40) :: fname

   INTEGER :: dX, dY, t, X0, X1, Y0, Y1
   INTEGER :: st
   
!  MPI variables
   CHARACTER(MPI_MAX_PROCESSOR_NAME) :: server
   INTEGER :: ierr, my_rank, numprocs, ls, dest, tag
   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
      CALL CPU_TIME(tempo1)

      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
!         WRITE(*,*) "Process ",my_rank," ALLOCATE(iters(",0,":",XYdots-1,",",&
!    &           0,":",XYdots-1,")"

      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
      call MPI_BCAST(Sreal,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(Simag,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(XYdots,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(Range,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)
      call MPI_BCAST(NITER,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)      
   WRITE(*,*) "Running MPI Mandel on server ",TRIM(server),". Input data are: ", &
   &           Sreal, Simag, XYdots, Range, procs, NITER

   ! Every process executes
   procs = numprocs   ! procs data required for compatibility reasons

   ! ---> Compute complex spacing between pixels <--- 
   XINC = Range / real(XYdots)  
   YINC = Range / real(XYdots)  
  
!  Let's decompose on 2nd dimension 
!   dX = ( XYdots + (Procs - 1) ) / Procs
   dY = ( XYdots + (Procs - 1) ) / Procs
!!   DO t = 0, Procs-1
!       X0 = dX * t
!       X1 = dX * ( t + 1)
!       if ( X1 > XYdots ) X1 = XYdots
       
       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,Y0:Y1-1),STAT=st)
         ALLOCATE(iters(0:XYdots-1,0:XYdots-1),STAT=st)
!!          WRITE(*,*) "Process ",my_rank," ALLOCATE(iters(",0,":",XYdots-1,",",&
!!          &           Y0,":",Y1-1,")"
!          IF ( st /= 0 ) THEN
!   	         WRITE(*,*) "Error allocating ITERS(",XYdots,",",XYdots,")"
!   	         LetsGo = .FALSE.
!   	         CALL Terminate(my_rank)
!          END IF
       END IF
       CALL calcola(X0, X1, Y0, Y1)
       IF ( my_rank /= 0 ) THEN
       !   Send data to master
          dest = 0; tag = 100
!!       PRINT*,"Process ",my_rank," sends iters(",X0,":", X1-1,",", Y0,":", Y1-1,")", &
!!  &            " to process ",dest
          call MPI_Send(iters(0,Y0), (XYdots*(Y1-Y0)), MPI_INTEGER, dest, &
                       tag, MPI_COMM_WORLD, ierr)
!!          PRINT*,"Process ",my_rank," after MPI_Send ierr = ",ierr
       END IF
!!   END DO
	
   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
!!       PRINT*,"Process ",my_rank," receives iters(",X0,":", X1-1,",", Y0,":", Y1-1,")", &
!!  &            " from process ",t
       tag = 100
       call MPI_Recv(iters(0,Y0), (XYdots*(Y1-Y0)), MPI_INTEGER, t, &
                       tag, MPI_COMM_WORLD, status, ierr)       
!!          PRINT*,"Process ",my_rank," after MPI_Recv ierr = ",ierr
      END DO
   
   CALL CPU_TIME(tempo2)
   
   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 :: pos
   INTEGER ::       IX      =  0   !  General purpose counter.     
   INTEGER ::       IY      =  0   !  General purpose counter.     
   INTEGER ::       IZ      =  0   !  General purpose counter.     
   REAL(8) ::     CA        = 0.0  !  Current coordinate. Real.     
   REAL(8) ::     CB        = 0.0  !  Current coordinate. Imaginary.    
   REAL(8) ::     ZA        = 0.0  !  Current itteration. Real.     
   REAL(8) ::     ZB        = 0.0  !  Current itteration. Imaginary.    
   REAL(8) ::     ZSIZE     = 0.0  !  Real distance of Z(A,B) from 0,0  
   REAL(8) ::     ZNEWA     = 0.0  !  Next itteration. Real.     
   REAL(8) ::     ZNEWB     = 0.0  !  Next itteration. Immaginary.     

   WRITE(*,*) "SHAPE(iters) = ",SHAPE(iters)
   WRITE(*,*) "LBOUND(iters) = ",LBOUND(iters)
   WRITE(*,*) "UBOUND(iters) = ",UBOUND(iters)
   
   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