program main

!*****************************************************************************80
!
!! stats_test() tests stats().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    08 May 2021
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer i
  real ( kind = rk ) inf
  integer n
  integer n2
  real ( kind = rk ) r8vec_max
  real ( kind = rk ) r8vec_min
  real ( kind = rk ) r8vec_sum
  real ( kind = rk ) x2(10)
  real ( kind = rk ) x_max
  real ( kind = rk ) x_max2
  real ( kind = rk ) x_mean
  real ( kind = rk ) x_mean2
  real ( kind = rk ) x_min
  real ( kind = rk ) x_min2
  real ( kind = rk ) x_std
  real ( kind = rk ) x_std2
  real ( kind = rk ) x_sum
  real ( kind = rk ) x_sum2
  real ( kind = rk ) x_var
  real ( kind = rk ) x_var2

  interface
    subroutine stats ( x, n, x_sum, x_min, x_mean, x_max, x_var, x_std )
      integer, parameter :: rk = kind ( 1.0D+00 )
      integer, optional :: n
      real ( kind = rk ), optional :: x
      real ( kind = rk ), optional :: x_max
      real ( kind = rk ), optional :: x_mean
      real ( kind = rk ), optional :: x_min
      real ( kind = rk ), optional :: x_std
      real ( kind = rk ), optional :: x_sum
      real ( kind = rk ), optional :: x_var
    end subroutine
  end interface

  call timestamp ( )
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'stats_test():'
  write ( *, '(a)' ) '  FORTRAN90 version'
  write ( *, '(a)' ) '  Test stats(), with the interface:'
  write ( *, '(a)' ) '    [n,sum,min,mean,max,var,std] = stats(x)'
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) '    Call 10 times with random values'
  write ( *, '(a)' ) '    Compare results with MATLAB vector calculation.'

  n2 = 10
  call random_number ( harvest = x2(1:n2) )

  call stats ( )
  do i = 1, n2
    write ( *, '(g14.6)' ) x2(i)
    call stats ( x = x2(i) )
  end do

  inf = huge ( x2(1) )
  call stats ( x = inf, n = n, x_sum = x_sum, x_min = x_min, &
    x_mean = x_mean, x_max = x_max, x_var = x_var, x_std = x_std )

  x_sum2 = r8vec_sum ( n2, x2 )
  x_min2 = r8vec_min ( n2, x2 )
  call r8vec_mean ( n2, x2, x_mean2 )
  x_max2 = r8vec_max ( n2, x2 )
  call r8vec_variance_sample ( n2, x2, x_var2 )
  call r8vec_std_sample ( n2, x2, x_std2 )

  write ( *, '(a)' ) ''
  write ( *, '(a,i4,a,i4)' ) '      n1    = ', n,      ', n2    = ', n2
  write ( *, '(a,g14.6,a,g14.6)' ) '      sum1  = ', x_sum,  ', sum2  = ', x_sum2
  write ( *, '(a,g14.6,a,g14.6)' ) '      min1  = ', x_min,  ', min2  = ', x_min2
  write ( *, '(a,g14.6,a,g14.6)' ) '      mean1 = ', x_mean, ', mean2 = ', x_mean2
  write ( *, '(a,g14.6,a,g14.6)' ) '      max1  = ', x_max,  ', max2  = ', x_max2
  write ( *, '(a,g14.6,a,g14.6)' ) '      var1  = ', x_var,  ', var2  = ', x_var2
  write ( *, '(a,g14.6,a,g14.6)' ) '      std1  = ', x_std,  ', std2  = ', x_std2
!
!  Terminate.
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'stats_test():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ''
  call timestamp ( )

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 August 2021
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 8 ) ampm
  integer d
  integer h
  integer m
  integer mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer n
  integer s
  integer values(8)
  integer y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2.2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end