!
! *****************************************************************
! *                                                               *
! *    Copyright (c) Digital Equipment Corporation, 1995          *
! *                                                               *
! *   All Rights Reserved.  Unpublished rights  reserved  under   *
! *   the copyright laws of the United States.                    *
! *                                                               *
! *   The software contained on this media  is  proprietary  to   *
! *   and  embodies  the  confidential  technology  of  Digital   *
! *   Equipment Corporation.  Possession, use,  duplication  or   *
! *   dissemination of the software and media is authorized only  *
! *   pursuant to a valid written license from Digital Equipment  *
! *   Corporation.                                                *
! *                                                               *
! *   RESTRICTED RIGHTS LEGEND   Use, duplication, or disclosure  *
! *   by the U.S. Government is subject to restrictions  as  set  *
! *   forth in Subparagraph (c)(1)(ii)  of  DFARS  252.227-7013,  *
! *   or  in  FAR 52.227-19, as applicable.                       *
! *                                                               *
! *****************************************************************
!
!
! HISTORY
!
!+
! HPF Mandelbrot program
! Written in DEC Fortran 90 by Digital Equipment Corporation
!
! Copyright (c) 1994 by Digital Equipment Corporation
! Author: Jonathan Harris
! Contributions from David LaFrance-Linden, David Loveman
!
! Compile this program as follows:
!   f90 -o mandelbrot mandelbrot.f90 simpleX.c -wsf -lX11 -lm -ldnet_stub
!
! Run the program by typing:
!   mandelbrot -peers N
!                     ^Replace N by the number of nodes in your farm, eg.
!                         mandelbrot -peers 8
!
!-

!+
! Timer routines
!-
module timer
contains

  !+
  ! cputime - Returns elapsed time in seconds since previous invocation;
  !           returns 0.0 the first time it is called.
  !-
  real function cputime()
    integer :: clock_count, clock_rate, clock_max
    integer, save :: last_count
    logical, save :: system_clock_not_yet_called = .true.

    call system_clock(count=clock_count,     &
                      count_rate=clock_rate,   &
                      count_max=clock_max)
    if (system_clock_not_yet_called) then
      system_clock_not_yet_called = .false.
      cputime = 0.0
    else if (clock_count > last_count) then
      cputime = real(clock_count-last_count)/real(clock_rate)
    else
      cputime = real(clock_max-last_count+clock_count)  &
                /real(clock_rate)
    end if
    last_count = clock_count
    return
  end function cputime

end module timer

!+
! Main program
!-
program mandelbrot
  use timer

  integer, parameter :: w=625, h=625           ! Width, height of image
  integer, parameter :: cx = w/2, cy = h/2     ! Parameters for computation
  double precision, parameter   :: radius = 2.0

  byte	             :: target(h,w)            ! Image being computed lives here
  !HPF$ distribute target(*,cyclic)

  double precision   :: ccr = 0.0D0, cci = 0.0D0 !Starting point for computation
  double precision   :: cstep = radius * 2 / min(w, h)
  integer            :: nmin=0, nmax=255       ! Min, max # of iters (colors)

  real               :: etime, dtime           ! Elapsed time

  logical            :: done                   ! Controls termination of loop

  interface
    ! Interface to mandel_val function
    pure byte function mandel_val(x, nmin, nmax)
      complex(KIND=KIND(0.0D0)), intent(in) ::   x
      integer, intent(in)  ::   nmin, nmax
    end function mandel_val

    ! Interface to X setup routine
    extrinsic (hpf_serial) subroutine MBSetupX(h, w, ccr, cci, cstep, nprocs)
      integer, intent(in)	:: h, w, nprocs
      double precision, intent(in)   :: ccr, cci, cstep
    end subroutine MBSetupX

    ! Interface to X display routine
    extrinsic (hpf_serial) subroutine MBdisplayX(target0)
      byte, intent(in)   :: target0(:,:)
    end subroutine MBdisplayX

    ! Interface to X input routine
    extrinsic (hpf_serial) subroutine MBinputX(ccr, cci, cstep, nmin, nmax, done)
      double precision, intent(inout)	:: ccr, cci, cstep
      integer, intent(inout)	:: nmin, nmax
      logical, intent(out)	:: done
    end subroutine MBinputX

  end interface

  ! Setup X window to display Mandelbrot picture
  call MBSetupX(h, w, ccr, cci, cstep, number_of_processors())

  !
  ! Main loop: Compute Mandelbrot image, get coordinates of next image, loop
  !
  done = .false.

  do while (.not. done)
  print "(/, 'Computing ',i4,' by ',i4,' Mandelbrot set at [',f6.3,', ',f6.3," &
        // "', ',e9.3,'] on ',i2,' processor Alpha farm')", &
        w, h, ccr, cci, cstep, number_of_processors()

  ! Compute Mandelbrot image: Call mandel_val function for each point in image

  etime=cputime()                               ! Start timer
  forall(ix = 1:w, iy = 1:h)                                  &
    target(ix,iy) = mandel_val(CMPLX(ccr + ((ix-cx)*cstep),   &
                                     cci + ((iy-cy)*cstep),   &
      	                             KIND=KIND(0.0D0)),       &
                               nmin, nmax)
  etime=cputime()
  print "('Compute time:		', f6.1, ' seconds')", etime

  ! Display Mandelbrot image we just computed
  dtime = cputime()
  call MBdisplayX(target)

  dtime = cputime()
  print "('Display time:		', f6.1, ' seconds')", dtime

  print "('Computation rate:	', f6.1, ' Mflops')", &
	real(9*sum(zext(target)))/(etime*1.0e6)

  ! Get coordinates of next image to compute
  call MBinputX(ccr, cci, cstep, nmin, nmax, done)

  end do

end program mandelbrot

!+
! mandel_val - Computes Mandelbrot function at a point
!-
pure byte function mandel_val(x, nmin, nmax)
  complex(KIND=KIND(0.0D0)), intent(in) ::   x
  integer, intent(in)                   :: nmin, nmax
	
  integer	                        :: n

!  The following commented-out code is replaced by the code below it for speed.
!
!  complex(KIND=KIND(0.0D0))  :: xtemp
!
!  n = 0
!  xtemp = x
!  do while (ABS(xtemp) < 2.0 .AND. (n < nmax))
!    xtemp = xtemp * xtemp + x
!    n = n + 1
!  end do

! Beginning of speedup

  real(kind=KIND(0.0D0)) :: xorgr, xorgi, xr, xi, xr2, xi2, rad2
  logical                :: keepgoing

  n = -1
  xorgr = REAL(x)
  xorgi = AIMAG(x)
  xr = xorgr
  xi = xorgi

  do
    n = n + 1
    xr2 = xr*xr
    xi2 = xi*xi
    xi = 2*(xr*xi) + xorgi
    keepgoing = n < nmax
    rad2 = xr2 + xi2
    xr = xr2 - xi2 + xorgr
    if (keepgoing .AND. (rad2 <= 4.0)) cycle
    exit
  end do

!  End of speedup version of code
	
  if (n >= nmax) then
    mandel_val = nmax-nmin
  else
    mandel_val = MOD(n, nmax-nmin)
  end if
	
end function mandel_val

