! Calculate the area of the Mandelbrot set by simple point counting.
! Because it is symmetric about the X-axis, we save time by doing
! that only for the upper half.  It offsets the points to correct for
! the double counting of the points on the axis.

! This uses process 0 as the root process, for reading from stdin;
! that usually works, but there is a later lecture on I/O.  

! Look for UseMPI to see the MPI extensions to the serial code.

PROGRAM Mandelbrot
    USE double
    USE mpi
    IMPLICIT NONE
    LOGICAL, PARAMETER :: UseMPI = .True.
    INTEGER, PARAMETER :: root = 0, tag_1 = 123, tag_2 = 456
    INTEGER :: maxiters, npoints, error, nprocs, myrank, count, rate
    LOGICAL :: finished = .False.
    REAL(KIND=DP) :: stamp, buffer_1(2), step, total = 0.0_DP
    COMPLEX(KIND=DP) :: buffer_2(2)

! Initialise MPI and find out our environment. 

    IF (UseMPI) THEN
        CALL MPI_Init(error)
        CALL MPI_Comm_size(MPI_COMM_WORLD,nprocs,error)
        CALL MPI_Comm_rank(MPI_COMM_WORLD,myrank,error)
        stamp = MPI_Wtime()
    ELSE
        nprocs = 1
        myrank = root
        CALL SYSTEM_CLOCK(count,rate)
        stamp = count/(rate+0.0D0)
    END IF

! Read in the parameters.

    IF (myrank == root) THEN
        READ *, maxiters, step, npoints
        IF (maxiters < 10) THEN
            PRINT *, 'Invalid value of MAXITERS'
            CALL MPI_Abort(MPI_COMM_WORLD,1,error)
        END IF
        IF (step < 10.0_DP*EPSILON(step) .OR. step > 0.1_DP) THEN
            PRINT *, 'Invalid value of STEP'
            CALL MPI_Abort(MPI_COMM_WORLD,1,error)
        END IF
        IF (npoints < 1 .OR. npoints > 4.0/((nprocs-1)*step)) THEN
            PRINT *, 'Invalid value of NPOINTS'
            CALL MPI_Abort(MPI_COMM_WORLD,1,error)
        END IF
        buffer_1(1) = maxiters
        buffer_1(2) = step
    END IF

! Distribute the parameters using broadcast, and call one procedure for
! the root process and another for all of the others.

    IF (UseMPI .AND. nprocs > 1) THEN
        CALL MPI_Bcast(buffer_1,2,MPI_DOUBLE_PRECISION,    &
            root,MPI_COMM_WORLD,error)
        maxiters = buffer_1(1)
        step = buffer_1(2)
        IF (myrank == root) THEN
            CALL Rootproc(total)
        ELSE
            CALL Otherproc
        END IF

! Do the job in serial mode.

    ELSE
        total = 0.0_DP
        DO
            CALL Iterate (finished, buffer_2)
            IF (finished) EXIT
            total = total + Shell(buffer_2(1),buffer_2(2))
        END DO
    END IF

! Print the result and finish.

    IF (myrank == root) THEN
        PRINT '(A,F0.3)', 'Area of Mandelbrot set is about ',    &
            2.0_DP*total
        IF (UseMPI) THEN
            PRINT '(A,F0.1," seconds")', 'Time taken was about ',    &
                MPI_Wtime()-stamp
        ELSE
            CALL SYSTEM_CLOCK(count,rate)
            PRINT '(A,F0.1," seconds")', 'Time taken was about ',    &
                count/(rate+0.0D0)-stamp
        END IF
    END IF
    IF (UseMPI) THEN
        CALL MPI_Finalize(error)
    END IF

CONTAINS



! This code is executed in the root process.  Note that it starts off
! by seeding all of the other processes, finishes by sending all of
! them an EOT indicator, and tops them up as they finish in between.
! There is little point in making the sends non-blocking for a small
! number of processes, though it might help for a large one.
!
! Note that it assumes external checking - i.e. that it is run on the
! root and that there is at least once domain per processor.

    SUBROUTINE Rootproc (total)
        REAL(KIND=DP), INTENT(OUT) :: total
        INTEGER :: status(MPI_STATUS_SIZE)
        INTEGER :: index, sent = 0, received = 0, error, n
        INTEGER, ALLOCATABLE :: requests(:)
        LOGICAL :: finished = .False.
        COMPLEX(KIND=DP) :: buffer_2(2)
        REAL(KIND=DP), ALLOCATABLE :: buffer_3(:)

        ALLOCATE(requests(0:nprocs-1),buffer_3(0:nprocs-1))
        total = 0.0_DP
        DO n = 0,nprocs-1
            IF (n /= root) THEN
                CALL Iterate(finished,buffer_2)
                CALL MPI_Send(buffer_2,2,MPI_DOUBLE_COMPLEX,    &
                    n,tag_1,MPI_COMM_WORLD,error)
                CALL MPI_Irecv(buffer_3(n),1,MPI_DOUBLE_PRECISION,    &
                    n,tag_2,MPI_COMM_WORLD,requests(n),error)
                sent = sent+1
            ELSE
                requests(n) = MPI_REQUEST_NULL
            END IF
        END DO
        DO WHILE (sent > received)
            CALL MPI_Waitany(nprocs,requests,index,status,error)
            received = received+1
            index = index-1
            IF (status(MPI_SOURCE) /= index) THEN
                PRINT *, 'Logic error in program'
                CALL MPI_Abort(MPI_COMM_WORLD,1,error)
            END IF
            total = total + buffer_3(index)
            IF (.NOT. finished) THEN
                CALL Iterate(finished,buffer_2)
                IF (finished) CYCLE
                CALL MPI_Send(buffer_2,2,MPI_DOUBLE_COMPLEX,    &
                    index,tag_1,MPI_COMM_WORLD,error)
                CALL MPI_Irecv(buffer_3(index),1,MPI_DOUBLE_PRECISION,    &
                    index,tag_2,MPI_COMM_WORLD,requests(index),error)
                sent = sent+1
            END IF
        END DO
        buffer_2 = CMPLX(0.0_DP,-1.0_DP,KIND=DP)
        DO n = 0,nprocs-1
            CALL MPI_Send(buffer_2,2,MPI_DOUBLE_COMPLEX,    &
                n,tag_1,MPI_COMM_WORLD,error)
        END DO
        DEALLOCATE(requests,buffer_3)
    END SUBROUTINE Rootproc



! This code is executed in each of the worker processes.

    SUBROUTINE Otherproc
        INTEGER :: status(MPI_STATUS_SIZE)
        INTEGER :: error
        COMPLEX(KIND=DP) :: buffer_2(2)
        REAL(KIND=DP) :: result

        DO
            CALL MPI_Recv(buffer_2,2,MPI_DOUBLE_COMPLEX,    &
               root,tag_1,MPI_COMM_WORLD,status,error)
            IF (AIMAG(buffer_2(1)) >= 0.0_DP) THEN
                result = Shell(buffer_2(1),buffer_2(2))
                CALL MPI_Send(result,1,MPI_DOUBLE_PRECISION,    &
                   root,tag_2,MPI_COMM_WORLD,error)
            ELSE
                EXIT
            END IF
        END DO
    END SUBROUTINE Otherproc



! Return a sequence of domains that cover the area we are integrating
! over.  Note that it returns the limits as an array of length 2.

    SUBROUTINE Iterate (finished, domain)
        IMPLICIT NONE
        LOGICAL, INTENT(OUT) :: finished
        COMPLEX(KIND=DP), INTENT(OUT) :: domain(2)
        COMPLEX(KIND=DP), SAVE :: initial = (-2.0_DP,0.0_DP)

        finished = (AIMAG(initial) >= 2.0_DP)
        IF (finished) RETURN
        domain = initial
        domain(2) = domain(2)+CMPLX(step*npoints,step*npoints,KIND=DP)
        IF (REAL(initial) < 2.0_DP) THEN
            initial = initial+step*npoints
        ELSE
            initial = CMPLX(-2.0_DP,AIMAG(initial)+step*npoints,KIND=DP)
        END IF
    END SUBROUTINE Iterate



! Scan over the domain we have been given, and return the area in
! the set.  Note that it is passed the limits as separate arguments.

    PURE FUNCTION Shell (lower, upper)
        IMPLICIT NONE
        REAL(KIND=DP) :: Shell
        COMPLEX(KIND=DP), INTENT(IN) :: lower, upper
        COMPLEX(KIND=DP) :: work

        Shell = 0.0_DP
        work = CMPLX(REAL(lower)+step/2.0_DP,    &
            AIMAG(lower)+step/2.0_DP,KIND=DP)
        DO WHILE (AIMAG(work) < AIMAG(upper))
            DO WHILE (REAL(work) < REAL(upper))
                IF (Kernel(work)) Shell = Shell+step**2
                work = work+step
            END DO
            work = CMPLX(REAL(lower)+step/2.0_DP,    &
                AIMAG(work)+step,KIND=DP)
        END DO
    END FUNCTION Shell



! See if this point is within the Mandelbrot set.

    PURE FUNCTION Kernel (value)
        IMPLICIT NONE
        LOGICAL :: Kernel
        COMPLEX(KIND=DP), INTENT(IN) :: value
        COMPLEX(KIND=DP) :: work
        INTEGER :: n

        work = value
        DO n = 1, maxiters
            work = work**2 + value
            IF (ABS(REAL(work)) > 2.0 .OR. ABS(AIMAG(work)) > 2.0) EXIT
        END DO
        Kernel = (ABS(WORK) <= 2.0) 
    END FUNCTION Kernel

END PROGRAM Mandelbrot
