! o separated type definitions out of arrays_module and into new
!   types_module
! o added loop1b, loop1c, loop1d
! o replaced explicit (1:nc) with whole array operations
!  > cleaner, no effect on timings
! o replaced (1,1:nc) with (1,:) in loop6 and loop7
!  > cleaner, no effect on timings
! o augmented tests of (2,nc) and (9,nc) arrays
! o augmented loop2 test

module JT_Clock_module
  implicit none

  private
  public :: JT_Clock

contains

  subroutine JT_Clock (time)

    ! Out:
    real, intent(out) :: time

    ! Local:
    real :: user, syst, cputot, ETIME

    cputot = ETIME(user, syst)
    time = user + syst

    return
  end subroutine JT_Clock

end module JT_Clock_module

module nc_module
  implicit none

  integer, parameter :: nc=3000
end module nc_module

module types_module
  use nc_module

  implicit none

  type scalar
     real :: x
  end type scalar

  type vector
     real, dimension(nc) :: x
  end type vector

  type two_scalars
     real :: x, y
  end type two_scalars

  type real_pointer
     real, pointer :: x
  end type real_pointer

  type nine_scalars
     real :: x1, x2, x3, x4, x5, x6, x7, x8, x9
  end type nine_scalars
end module types_module

module arrays_module
  use nc_module
  use types_module

  real, target ::  a(nc), b(nc), c(nc)
  real, allocatable ::  aa(:), ab(:), ac(:)

  integer :: ivec(nc)
  real :: a2(2,nc), c2(2,nc)
  real :: a9(9,nc), b9(9,nc), c9(9,nc)
  real, pointer :: pa(:), pb(:), pc(:), paa(:), pab(:), pac(:)
  real :: b2(2,nc)

  type(scalar) :: xa1(nc), xb1(nc), xc1(nc)
  type(vector) :: xa2, xb2, xc2
  type(two_scalars) :: xa3(nc), xb3(nc), xc3(nc)
  type(real_pointer) :: xa4(nc), xb4(nc), xc4(nc), sa4(nc), sb4(nc), sc4(nc)
  type(nine_scalars) :: xa9(nc), xb9(nc), xc9(nc)

end module arrays_module

subroutine loop1
  use arrays_module
  implicit none

  c(1:nc) = a(1:nc)*b(1:nc)

  return
end subroutine loop1

subroutine loop1a
  use arrays_module
  implicit none

  c = a*b

  return
end subroutine loop1a

subroutine loop1b (a, b, c, nc)
  implicit none
  integer :: nc
  real, dimension(nc) :: a, b, c

  c = a*b

  return
end subroutine loop1b

module loop1c_module
  implicit none

contains

  subroutine loop1c (a, b, c)
    real, dimension(:) :: a, b, c

    c = a*b

    return
  end subroutine loop1c
end module loop1c_module

module loop1d_module
  implicit none

contains

  subroutine loop1d (a, b, c)
    real, dimension(:), intent(in) :: a, b
    real, dimension(:), intent(out) :: c

    c = a*b

    return
  end subroutine loop1d
end module loop1d_module

subroutine loop2
  use arrays_module
  implicit none

  xc1%x = xa1%x*xb1%x

  return
end subroutine loop2

module loop2a_module
  implicit none

contains

  subroutine loop2a (a, b, c)
    use types_module
    type(scalar), dimension(:), intent(in) :: a, b
    type(scalar), dimension(:), intent(out) :: c

    c%x = a%x*b%x

    return
  end subroutine loop2a
end module loop2a_module

subroutine loop3
  use arrays_module
  implicit none

  xc2%x = xa2%x*xb2%x

  return
end subroutine loop3

subroutine loop4
  use arrays_module
  implicit none

  xc3%x = xa3%x*xb3%x

  return
end subroutine loop4

subroutine loop5
  use arrays_module
  implicit none
  integer :: j

  do  j=1,nc
     xc4(j)%x = xa4(j)%x*xb4(j)%x
  end do

  return
end subroutine loop5

subroutine loop6
  use arrays_module
  implicit none

  c2(2,:) = a2(2,:)*b2(2,:)

  return
end subroutine loop6

module loop6a_module
  implicit none

contains

  subroutine loop6a (a, b, c)
    real, dimension(:,:), intent(in) :: a, b
    real, dimension(:,:), intent(out) :: c

    c(2,:) = a(2,:)*b(2,:)

    return
  end subroutine loop6a
end module loop6a_module

subroutine loop7
  use arrays_module
  implicit none

  c9(2,:) = a9(2,:)*b9(2,:)

  return
end subroutine loop7

subroutine loop8
  use arrays_module
  implicit none

  xc9%x2=xa9%x2*xb9%x2

  return
end subroutine loop8

subroutine loop9
  use arrays_module
  implicit none

  pc = pa*pb

  return
end subroutine loop9

subroutine loop10
  use arrays_module
  implicit none

  c = a(ivec)*b(ivec)

  return
end subroutine loop10

subroutine loop11
  use arrays_module
  implicit none
  integer :: j

  do j = 1,nc
     sc4(j)%x = sa4(j)%x*sb4(j)%x
  end do

  return
end subroutine loop11

subroutine loop12
  use arrays_module
  implicit none

  ac = aa*ab

  return
end subroutine loop12

subroutine loop13 (a, b, c, nc)
  implicit none
  integer, intent(in) :: nc
  real, intent(in) :: a(nc), b(nc)
  real, intent(out) :: c(nc)

  c = a*b

  return
end subroutine loop13

module loop14_module
  implicit none

contains

  subroutine loop14 (a, b, c, nc)
    integer, intent(in) :: nc
    real, intent(in) ::  a(:), b(:)
    real, intent(out) :: c(:)

    c = a*b
  
    return
  end subroutine loop14
end module loop14_module

subroutine loop15
  use arrays_module
  implicit none

  pac = paa*pab

  return
end subroutine loop15

module loop16_module
  implicit none

contains

  subroutine loop16 (pa, pb, pc, nc)
    integer, intent(in) :: nc
    real, pointer :: pa(:), pb(:), pc(:)

    pc = pa*pb
  
    return
  end subroutine loop16
end module loop16_module

subroutine empty (a, b, c, nc)
  implicit none
  integer :: nc
  real :: a(nc), b(nc), c(nc)

  ! c(1:nc) = a(1:nc)*b(1:nc)
end subroutine empty

program speed
  use nc_module
  use arrays_module
  use JT_Clock_module
  use loop1c_module
  use loop1d_module
  use loop2a_module
  use loop6a_module
  use loop14_module
  use loop16_module

  implicit none

  integer, parameter :: nloop=3000
  integer :: i, iii, j, iswap
  real :: rannum, start, end, time

  ! Program to time several different array calculations

  a = 0.7
  a2 = 0.7
  b = 0.8
  b2 = 0.8

  ! Load the derived type arrays and assign pointers
  xa1%x = a
  xa2%x = a
  xa3%x = a
  xa9%x1 = a
  xb1%x = b
  xb2%x = b
  xb3%x = b
  xb9%x1 = b

  pa => a
  pb => b
  pc => c

  do i=1,nc
     xa4(i)%x => a(i)
     xb4(i)%x => b(i)
     xc4(i)%x => c(i)
     ivec(i) = i
  end do

  ! Assign random order in the indirect indices.
  do i=1,nc
     call random_number(rannum)
     j = nint(real((nc-1)*rannum)) + 1
     iswap = ivec(i)
     ivec(i) = ivec(j)
     ivec(j) = iswap
  end do

  do i=1,nc
     sa4(i)%x => a(ivec(i))
     sb4(i)%x => b(ivec(i))
     sc4(i)%x => c(i)
  end do

  allocate (aa(nc), ab(nc), ac(nc))
  allocate (paa(nc), pab(nc), pac(nc))
  aa = a
  ab = b

  ! The product of the array dimension "nc" and the number
  ! of calls to each subroutine "nloop" should be large enough to
  ! provide good statistics. "nc" is itself large enough that 
  ! the relative cost of calling the subroutine is insignificant.

  do iii=1,1

     print *
     print *, 'Results of element by element multiplication of 2 arrays'
     print *

     call JT_Clock (start)
     do i = 1,nloop
        call loop1
     end do
     call JT_Clock (end)
     print *, 'Normal arrays via module, explicit indexing = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1a
     end do
     call JT_Clock (end)
     print *, 'Normal arrays via module, whole array operation = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1b (a, b, c, nc)
     end do
     call JT_Clock (end)
     print *, 'Normal arrays via arg list = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1c (a, b, c)
     end do
     call JT_Clock (end)
     print *, 'Normal arrays via arg list, explicit interface = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1d (a, b, c)
     end do
     call JT_Clock (end)
     print *, 'Normal arrays via arg list, explicit interface, intent = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop6
     end do
     call JT_Clock (end)
     print *, ' (2,nc) arrays = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1d (a2(2,:), b2(2,:), c2(2,:))
     end do
     call JT_Clock (end)
     print *, ' (2,nc) arrays, segment passed in = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop6a (a2, b2, c2)
     end do
     call JT_Clock (end)
     print *, ' (2,nc) arrays, arrays passed in = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop7
     end do
     call JT_Clock (end)
     print *, ' (9,nc) arrays = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop1d (a9(2,:), b9(2,:), c9(2,:))
     end do
     call JT_Clock (end)
     print *, ' (9,nc) arrays, segment passed in = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop6a (a9, b9, c9)
     end do
     call JT_Clock (end)
     print *, ' (9,nc) arrays, arrays passed in = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop2
     end do
     call JT_Clock (end)
     print *, 'Array of derived type with 1 real component = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop2a (xa1, xb1, xc1)
     end do
     call JT_Clock (end)
     print *, 'Array of derived type with 1 real component, passed in = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop4
     end do
     call JT_Clock (end)
     print *, 'Array of derived type with 2 real components = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop8
     end do
     call JT_Clock (end)
     print *, 'Array of derived type with 9 real components = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop3
     end do
     call JT_Clock (end)
     print *, 'Derived type containing xr(nc) = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop5
     end do
     call JT_Clock (end)
     print *, 'Array of derived type with 1 pointer = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop9
     end do
     call JT_Clock (end)
     print *, 'Normal pointer to basic arrays = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop10
     end do
     call JT_Clock (end)
     print *, 'Indirect addressing = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop11
     end do
     call JT_Clock (end)
     print *, 'Above with array of pointers = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop12
     end do
     call JT_Clock (end)
     print *, 'Allocated arrays = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop15
     end do
     call JT_Clock (end)
     print *, 'Allocated pointers = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop13 (a, b, c, nc)
     end do
     call JT_Clock (end)
     print *, 'Arrays in argument list = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop13 (pa, pb, pc, nc)
     end do
     call JT_Clock (end)
     print *, 'Pointers in argument list without explicit interface = ', end-start,' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop16 (pa, pb, pc, nc)
     end do
     call JT_Clock (end)
     print *, 'Pointers in argument list with explicit interface = ', end-start,' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop13 (aa, ab, ac, nc)
     end do
     call JT_Clock (end)
     print *, 'Allocated arrays in argument list = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call loop14 (a, b, c, nc)
     end do
     call JT_Clock (end)
     print *, 'Arrays in argument list with explicit interface = ', end-start, ' seconds'

     call JT_Clock (start)
     do i = 1,nloop
        call empty (a, b, c, nc)
     end do
     call JT_Clock (end)
     print *, 'Empty subroutine with above arguments = ', end-start, ' seconds'

     print *
     print *, '-----------------------------------------------'
     print *
  end do
  print *, c(nc), xc1(nc)%x, xc2%x(nc)
  print *, xc3(nc)%x, xc4(nc)%x
  stop

end program speed