! ++=======================================++
! ||                                       ||
! ||    THIS IS THE ALTERNATIVE VERSION    ||
! ||                                       ||
! ++=======================================++




! 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 version is a bit fancier, and shows what can be done by simple
! tuning.  It does that by attempting to balance the work between the
! processors.  Note that it is a memory hog compared to the other
! version - TANSTAAFL.

! 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
    INTEGER :: buffer_1(2), maxiters, error, nprocs, myrank,    &
        limit, index, count, rate, i, n
    REAL(KIND=DP) :: stamp, step, x
    COMPLEX(KIND=DP), ALLOCATABLE :: buffer_2(:), buffer_3(:)
    COMPLEX(KIND=DP) :: work

! 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
        stamp = count/(rate+0.0D0)
    END IF

! Read in the parameters.  

    IF (myrank == root) THEN
        READ *, maxiters, step
        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
        buffer_1(1) = maxiters
    END IF

! Allocate the second transfer array to be the size we actually need
! for the scatter, rounded up to a multiple of nprocs.

    limit = ((2*CEILING(2.0_DP/step)**2-1)/nprocs+1)*nprocs
    IF (myrank == root) THEN
        ALLOCATE(buffer_2(limit))

! Fill it with the points within the unit circle.

        index = 0
        work = CMPLX(step/2.0_DP-2.0_DP,step/2.0_DP,KIND=DP)
        DO WHILE (AIMAG(work) < 2.0_DP)
            DO WHILE (REAL(work) < 2.0_DP)
                IF (REAL(work)**2+AIMAG(work)**2 < 4.0_DP) THEN
                    index = index+1
                    IF (index > limit) THEN
                        PRINT *, 'Oops - an internal error'
                        CALL MPI_Abort(MPI_COMM_WORLD,1,error)
                    END IF
                    buffer_2(index) = work
                END IF
                work = work+step
            END DO
            work = CMPLX(step/2.0_DP-2.0_DP,AIMAG(work)+step,KIND=DP)
        END DO

! Pad to a multiple of nprocs.

        DO WHILE (MOD(index,nprocs) /= 0)
            index = index+1
            buffer_2(index) = 2.0_DP
        END DO
        buffer_1(2) = index/nprocs

! Randomise its contents.

       DO i = 2,index
            CALL RANDOM_NUMBER(x)
            n = FLOOR(x*i+1)
            IF (n < i) THEN
                work = buffer_2(i)
                buffer_2(i) = buffer_2(n)
                buffer_2(n) = work
            END IF
        END DO
    ELSE
        ALLOCATE(buffer_2(1))    ! This is not actually used
    END IF

! Distribute the data, using broadcast for the parameters and scatter
! for the domains.  Note that the allocate the receiving buffer after
! receiving its size.

    IF (UseMPI) THEN
        CALL MPI_Bcast(buffer_1,2,MPI_INTEGER,    &
            root,MPI_COMM_WORLD,error)
        maxiters = buffer_1(1)
        index = buffer_1(2)
        ALLOCATE(buffer_3(index))
        CALL MPI_Scatter(buffer_2,index,MPI_DOUBLE_COMPLEX,    &
            buffer_3,index,MPI_DOUBLE_COMPLEX,root,MPI_COMM_WORLD,error)
    ELSE
        buffer_3 = buffer_2
    END IF

! Scan over the points we have been given and calculate the number
! of interior points.

    n = 0
    DO i = 1,index
        IF (Kernel(buffer_3(i))) n = n+1
    END DO
    buffer_1(1) = n

! Accumulate the total, print the result and finish.  Note that this
! buffer use is correct Fortran.

    IF (UseMPI) THEN
        CALL MPI_Reduce(buffer_1(1),buffer_1(2),1,MPI_INTEGER,    &
            MPI_SUM,root,MPI_COMM_WORLD,error)
    ELSE
        buffer_1(2) = buffer_1(1)
    END IF
    IF (myrank == root) THEN
        PRINT '(A,F6.3)', 'Area of Mandelbrot set is about',    &
            2.0_DP*buffer_1(2)*step**2
        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

! 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
