!
! *****************************************************************
! *                                                               *
! *    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
!
! $RCSfile: pi_example.f90,v $ $Revision: 1.1.2.6 $ (DEC) $Date: 1995/02/21 22:07:56 $
!

module timer_module

  integer, private :: starting_count
  logical, private :: system_clock_not_yet_called = .true.
  real, private    :: time_then = 0.0
  real, private    :: time_since_last_call
  private cputime

contains

  subroutine print_wallclock_time()
    character*8    :: date    ! ccyymmdd
    character*10   :: time    ! hhmmss.sss
    call date_and_time(date,time)
    print *, date(5:6)//'/'//date(7:8)//'/'//date(3:4)//   &
       '  '//time(1:2)//':'//time(3:4)//':'//time(5:10)
    return
  end subroutine print_wallclock_time

  real function cputime()
! This is a private function which calculates
! the number of seconds that the cpu devoted to the program

    integer :: current_count, clock_rate, maximum_count

    call system_clock(count=current_count,     &
                      count_rate=clock_rate,   &
                      count_max=maximum_count)
    if (system_clock_not_yet_called) then
      starting_count = current_count
      system_clock_not_yet_called = .false.
      cputime = 0.0
    else if (current_count > starting_count) then
      cputime = real(current_count-starting_count)/real(clock_rate)
    else
      cputime = real(maximum_count-starting_count+current_count)  &
                /real(clock_rate)
    end if
    time_since_last_call = cputime - time_then
    time_then = cputime
!    print *, '>>CPUTIME=', cputime
    return
  end function cputime

  subroutine begin_cpu_timer()
    real :: t
    t = cputime()
    return
  end subroutine begin_cpu_timer

  real function elapsed_cpu_time()
    real :: t
    t = cputime()
    elapsed_cpu_time = time_since_last_call
    return
  end function elapsed_cpu_time

end module timer_module

program pi_example
! compute the value of pi by numeric integration
  use timer_module
  implicit none
  !hpf$ nosequence

  integer, parameter   :: n = 100000     ! number of rectangles
  real, parameter      :: h = 1.0 / n   ! width of a rectangle
  integer, parameter   :: number_of_iterations = 15

  real                 :: time

  !hpf$ template t(n)
  !hpf$ distribute t(block)

  integer              :: i             ! forall index
  integer              :: j             ! do index
  real                 :: x             ! arg to arith stmt fcn

  double precision     :: pi
  double precision, dimension(n) :: rect_area
    !hpf$ align with t :: rect_area

  double precision ff
  ff(x) =  4 / (1.0 + x*x)              ! function to be integrated

  call print_wallclock_time()
  call begin_cpu_timer()

! first variation - arithmetic statement function 
  do j = 1, number_of_iterations
    forall (i=1:n) rect_area(i) = h * ff(h*(i-0.5))
    pi = sum(rect_area)
  end do

  time = elapsed_cpu_time()
  print *, ' '
  print *, 'value computed = ', pi
  print *, 'for arithmetic statement function called in forall'
  print 50 , '   elapsed time = ', time
50 format (a18, f12.8)
  call begin_cpu_timer()


! second variation - manually inlined arithmetic statement function
  do j = 1, number_of_iterations
    forall (i=1:n) 
      rect_area(i) = h * (4 / (1.0  +  h*(i-0.5) * h*(i-0.5)))
    end forall
    pi = sum(rect_area)
  end do

  time = elapsed_cpu_time()
  print *, ' '
  print *, 'value computed = ', pi
  print *, 'for arithmetic statement function inlined in forall'
  print 50, '   elapsed time = ', time
  call begin_cpu_timer()


! third variation - internal function
  do j = 1, number_of_iterations
    forall (i=1:n) rect_area(i) = h * f(h*(i-0.5))
    pi = sum(rect_area)
  end do

  time = elapsed_cpu_time()
  print *, ' '
  print *, 'value computed = ', pi
  print *, 'for internal pure function called in forall'
  print 50, '   elapsed time = ', time
  call begin_cpu_timer()


! fourth variation - scalar do loop with function inlined
  do j = 1, number_of_iterations
    pi = 0.0
    do i=1, n
      pi = pi + h * (4 / (1.0  +  h*(i-0.5) * h*(i-0.5)))
    end do
  end do

  time = elapsed_cpu_time()
  print *, ' '
  print *, 'value computed = ', pi
  print *, 'for scalar do loop with function inlined'
  print 50, '   elapsed time = ', time

  call print_wallclock_time()

contains

  pure double precision  function f(x)  ! function to be integrated
    real, intent(in) :: x
    f = 4 / (1.0 + x*x)
  end function f

end program pi_example
