program cluster_energy
!
!*******************************************************************************
!
!! CLUSTER_ENERGY computes the minimal cluster energy for a set of data.
!
!
!  Discussion:
!
!    The current code only generates the data internally, and does
!    not read an input file of data, which would be more useful.
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
  implicit none
!
  integer dim_num
!
  integer c_num
  real, allocatable, dimension ( :, : ) :: center
  integer c_hi
  integer c_lo
  real, allocatable, dimension ( : ) :: e
  integer ierror
  integer it_hi
  integer it_lo
  real, allocatable, dimension ( :, :) :: point
  integer point_dist
  integer point_num
  real, allocatable, dimension ( : ) :: r_max
  real, allocatable, dimension ( : ) :: r_min
  integer seed
  real t
!
  call timestamp ( )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CLUSTER_ENERGY'
  write ( *, '(a)' ) '  Compute the cluster energy of a set of points.'
  write ( *, '(a)' ) ' '

  seed = 1234567

  call random_initialize ( seed )

  call i_input ( 'Enter spatial dimension DIM_NUM', dim_num, ierror )

  allocate ( r_min(1:dim_num) )
  allocate ( r_max(1:dim_num) )

  call input_get ( point_dist, dim_num, r_min, r_max, point_num, c_lo, &
    c_hi, it_lo, it_hi )

  call input_check ( point_dist, dim_num, r_min, r_max, point_num, c_lo, &
    c_hi, it_lo, it_hi  )

  call input_print ( point_dist, dim_num, r_min, r_max, point_num, c_lo, &
    c_hi, it_lo, it_hi  )

  allocate ( point(1:dim_num,1:point_num) )

  call point_generate ( point_dist, dim_num, r_min, r_max, point_num, point )

  if ( .false. ) then
    call point_print ( dim_num, point_num, point )
  end if

  allocate ( e(c_lo:c_hi) )
  allocate ( center(1:dim_num,1:c_hi) )

  do c_num = c_lo, c_hi

    call cluster_iteration ( dim_num, point_num, c_num, point, r_min, r_max, &
      center, e(c_num), it_lo, it_hi  )

  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Energy table:'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Clusters  Energy  Energy/point'
  write ( *, '(a)' ) ' '

  do c_num = c_lo, c_hi
    t = e(c_num) / real ( point_num )
    write ( *, '(i4,3g12.6)' ) c_num, e(c_num), t, sqrt ( t )
  end do

  deallocate ( center )
  deallocate ( e )
  deallocate ( point )
  deallocate ( r_max )
  deallocate ( r_min )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CLUSTER_ENERGY'
  write ( *, '(a)' ) '  Normal end of execution.'

  stop
end
subroutine cluster_iteration ( dim_num, point_num, c_num, point, r_min, r_max, &
  center, e, it_lo, it_hi  )
!
!*******************************************************************************
!
!! CLUSTER_ITERATION seeks the minimal energy of a cluster of a given size.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer C_NUM, the number of clusters.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Output, real CENTER(DIM_NUM,C_NUM), the centers associated with
!    the minimal energy clustering.
!
!    Output, real E, the energy associated with the minimal energy clustering.
!
!    Input, integer IT_LO, IT_HI, the least and most number of energy
!    iterations.
!
  implicit none
!
  integer c_num
  integer dim_num
  integer point_num
!
  real, dimension ( dim_num, c_num ) :: center
  real, dimension ( dim_num, c_num ) :: centroid
  integer, dimension ( point_num ) :: cluster
  real e
  real e_new
  real e_new_new
  integer i
  integer it
  integer it_hi
  integer it_lo
  integer j
  integer missed
  real, dimension ( dim_num, point_num ) :: point
  real, dimension ( dim_num ) :: r
  real, dimension ( dim_num ) :: r_max
  real, dimension ( dim_num ) :: r_min
  integer, dimension ( c_num ) :: subs
!
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of clusters = ', c_num
!
!  Initialize the centers randomly.
!
  do j = 1, c_num
    call random_number ( harvest = r(1:dim_num) )
    center(1:dim_num,j) = ( 1.0E+00 - r(1:dim_num) ) * r_min(1:dim_num) &
                                    + r(1:dim_num)   * r_max(1:dim_num)
  end do
!
!  Initialize the clusters randomly.
!
  subs(1:c_num) = 0
  do i = 1, point_num
    call i_random ( 1, c_num, cluster(i) )
    j = cluster(i)
    subs(j) = subs(j) + 1
  end do

  call energy ( dim_num, point_num, c_num, point, center, cluster, e )

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  Initial energy =                        ', e

  do j = 1, c_num
    write ( *, '(2i5,6f10.4)' ) j, subs(j), center(1:dim_num,j)
  end do

  do it = 1, it_hi
!
!  #1:
!  Assign each point to the cluster of its nearest center.
!
    do i = 1, point_num
      call nearest_point ( dim_num, c_num, center, point(1,i), j )
      cluster(i) = j
    end do

    call energy ( dim_num, point_num, c_num, point, center, cluster, e )

    write ( *, '(a,g14.6)' ) '  Energy with old centers, new clusters = ', e
!
!  #2:
!  Determine the centroids of the clusters.
!
    centroid(1:dim_num,1:c_num) = 0.0E+00
    subs(1:c_num) = 0

    do i = 1, point_num
      j = cluster(i)
      subs(j) = subs(j) + 1
      centroid(1:dim_num,j) = centroid(1:dim_num,j) + point(1:dim_num,i)
    end do

    missed = 0

    do j = 1, c_num

      if ( subs(j) /= 0 ) then
        centroid(1:dim_num,j) = centroid(1:dim_num,j) / real ( subs(j) )
      else
        missed = missed + 1
        call random_number ( harvest = r(1:dim_num) )
        centroid(1:dim_num,j) = ( 1.0E+00 - r(1:dim_num) ) * r_min(1:dim_num) &
                                          + r(1:dim_num)   * r_max(1:dim_num)

      end if

    end do

    if ( missed /= 0 ) then
      write ( *, '(a,i6)' ) 'Number of empty clusters was ', missed
    end if

    call energy ( dim_num, point_num, c_num, point, centroid, cluster, e_new )

    write ( *, '(a,g14.6)' ) '  Energy with new centers, new clusters = ', e_new

    do j = 1, c_num
      write ( *, '(2i6,6f10.4)' ) j, subs(j), center(1:dim_num,j)
      write ( *, '(12x,6f10.4)' )              centroid(1:dim_num,j)
    end do
!
!  Update the centers.
!
    center(1:dim_num,1:c_num) = centroid(1:dim_num,1:c_num)

  end do

  return
end
subroutine energy ( dim_num, point_num, c_num, point, center, cluster, e )
!
!*******************************************************************************
!
!! ENERGY computes the energy of a given clustering.
!
!
!  Discussion:
!
!    The energy of the clustering is the sum of the energy of each cluster.
!
!    The sum of a cluster is the sum of the squares of the distances of 
!    each point in the cluster to the center point of the cluster.
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer C_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Input, real CENTER(DIM_NUM,C_NUM), the center points.
!
!    Input, integer CLUSTER(POINT_NUM), indicates the cluster to which
!    each data point belongs.
!
!    Output, real E, the energy of the clustering.
!
  implicit none
!
  integer c_num
  integer dim_num
  integer point_num
!
  real, dimension ( dim_num, c_num ) :: center
  integer, dimension ( point_num ) :: cluster
  real e
  integer i
  integer j
  real, dimension ( dim_num, point_num ) :: point
!
  e = 0.0E+00
  do i = 1, point_num
    j = cluster(i)
    e = e + sum ( ( center(1:dim_num,j) - point(1:dim_num,i) )**2 )
  end do

  return
end
subroutine i_input ( string, value, ierror )
!
!*******************************************************************************
!
!! I_INPUT prints a prompt string and reads an integer from the user.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Output, integer VALUE, the value input by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  integer ierror
  character ( len = * ) string
  integer value
!
  ierror = 0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )

  read ( *, *, iostat = ierror ) value

  return
end
subroutine i_random ( ilo, ihi, i )
!
!*******************************************************************************
!
!! I_RANDOM returns a random integer in a given range.
!
!
!  Modified:
!
!    23 September 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ILO, IHI, the minimum and maximum acceptable values.
!
!    Output, integer I, the randomly chosen integer.
!
  implicit none
!
  integer i
  integer ihi
  integer ilo
  real r
  real, parameter :: rhi = 1.0E+00
  real, parameter :: rlo = 0.0E+00
!
  call r_random ( rlo, rhi, r )

  i = ilo + int ( r * real ( ihi + 1 - ilo ) )
!
!  In case of oddball events at the boundary, enforce the limits.
!
  i = max ( i, ilo )
  i = min ( i, ihi )

  return
end
subroutine i_range_input ( string, value1, value2, ierror )
!
!*******************************************************************************
!
!! I_RANGE_INPUT reads a pair of integers from the user, representing a range.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Output, integer VALUE1, VALUE2, the values entered by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  integer ierror
  character ( len = * ) string
  integer value1
  integer value2
!
  ierror = 0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )
  read ( *, *, iostat = ierror ) value1, value2

  return
end
subroutine input_check ( point_dist, dim_num, r_min, r_max, point_num, &
  c_lo, c_hi, it_lo, it_hi )
!
!*******************************************************************************
!
!! INPUT_CHECK performs some simple checks on the problem data.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_DIST, the point distribution to use.
!    1, use a uniform random distribution.
!    2, use a uniform grid of points.  (This hasn't been set up properly
!       except for 1D!).
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer C_LO, C_HI, the minimum and maximum number of
!    clusters to consider.
!
!    Input, integer IT_LO, IT_HI, the least and most number of energy
!    iterations.
!
  implicit none
!
  integer dim_num
!
  integer c_hi
  integer c_lo
  integer it_hi
  integer it_lo
  integer point_dist
  integer point_num
  real r_max(dim_num)
  real r_min(dim_num)
!
  if ( point_dist < 1 .or. point_dist > 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) &
      '  Point distribution option less than 1 or greater than 2.'
    write ( *, '(a,i6)' ) '  POINT_DIST = ', point_dist
    stop
  end if

  if ( dim_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  Number of dimensions DIM_NUM < 1.'
    write ( *, '(a,i6)' ) '  DIM_NUM = ', dim_num
    stop
  end if

  if ( any ( r_min(1:dim_num) >= r_max(1:dim_num) ) ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  For at least one index I, R_MIN(I) >= R_MAX(I).'
    stop
  end if

  if ( point_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  Number of points POINT_NUM < 1.'
    write ( *, '(a,i6)' ) '  POINT_NUM = ', point_num
    stop
  end if

  if ( c_lo > c_hi ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  CLO_NUM > CHI_NUM.'
    write ( *, '(a,i6)' ) '  CLO_NUM = ', c_lo
    write ( *, '(a,i6)' ) '  CHI_NUM = ', c_hi
    stop
  end if

  if ( c_hi > point_num ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  CHI_NUM > POINT_NUM.'
    write ( *, '(a,i6)' ) '  CHI_NUM = ', c_hi
    write ( *, '(a,i6)' ) '  POINT_NUM = ', point_num
    stop
  end if

  if ( it_lo < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  IT_LO < 1.'
    write ( *, '(a,i6)' ) '  IT_LO = ', it_lo
    stop
  end if

  if ( it_lo > it_hi ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'INPUT_CHECK - Fatal error!'
    write ( *, '(a)' ) '  IT_LO > IT_HI.'
    write ( *, '(a,i6)' ) '  IT_LO = ', it_lo
    write ( *, '(a,i6)' ) '  IT_HI = ', it_hi
    stop
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INPUT_CHECK - The input data passed all checks.'

  return
end
subroutine input_get ( point_dist, dim_num, r_min, r_max, point_num, &
  c_lo, c_hi, it_lo, it_hi )
!
!*******************************************************************************
!
!! INPUT_GET gets the value of certain input quantities from the user.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer POINT_DIST, the point distribution to use.
!    1, use a uniform random distribution.
!    2, use a uniform grid of points.  (This hasn't been set up properly
!       except for 1D!).
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Output, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Output, integer POINT_NUM, the number of data points.
!
!    Output, integer C_LO, C_HI, the minimum and maximum number of
!    clusters to consider.
!
!    Output, integer IT_LO, IT_HI, the least and most number of energy
!    iterations.
!
  implicit none
!
  integer dim_num
!
  integer c_hi
  integer c_lo
  integer ierror
  integer it_hi
  integer it_lo
  integer point_dist
  integer point_num
  real, dimension ( dim_num ) :: r_max
  real, dimension ( dim_num ) :: r_min
!
  call rvec_range_input ( 'Enter lower and upper region bounds:', &
    dim_num, r_min, r_max, ierror )

  call i_input ( 'Enter number of data points POINT_NUM', point_num, ierror )

  call i_input ( 'Enter 1 for uniform random data, 2 for uniform grid', &
    point_dist, ierror )

  call i_range_input ( 'Enter lower and upper number of clusters', c_lo, &
    c_hi, ierror )

  call i_range_input ( 'Enter lower and upper number of energy iterations', &
    it_lo, it_hi, ierror )

  return
end
subroutine input_print ( point_dist, dim_num, r_min, r_max, point_num, &
  c_lo, c_hi, it_lo, it_hi )
!
!*******************************************************************************
!
!! INPUT_PRINT prints the value of the input quantities.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_DIST, the point distribution to use.
!    1, use a uniform random distribution.
!    2, use a uniform grid of points.  (This hasn't been set up properly
!       except for 1D!).
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer C_LO, C_HI, the minimum and maximum number of
!    clusters to consider.
!
!    Input, integer IT_LO, IT_HI, the least and most number of energy
!    iterations.
!
  implicit none
!
  integer dim_num
!
  integer c_hi
  integer c_lo
  integer i
  integer it_hi
  integer it_lo
  integer point_dist
  integer point_num
  real r_max(dim_num)
  real r_min(dim_num)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'INPUT_PRINT:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of spatial dimensions = ', dim_num
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '  Lower and upper region bounds:'
  write ( *, '(a)' ) ' '
  do i = 1, dim_num
    write ( *, '(4x,2f10.6)' ) r_min(i), r_max(i)
  end do

  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of data points =              ', point_num
  write ( *, '(a,i6)' ) '  Lowest number of clusters =  	', c_lo
  write ( *, '(a,i6)' ) '  Highest number of clusters =         ', c_hi
  write ( *, '(a,i6)' ) '  Least number of energy iterations =  ', it_lo
  write ( *, '(a,i6)' ) '  Most number of energy iterations =   ', it_hi
  write ( *, '(a)' ) ' '
  if ( point_dist == 1 ) then
    write ( *, '(a)' ) '  Data points will be uniform random values.'
  else if ( point_dist == 2 ) then
    write ( *, '(a)' ) '  Data points will be a uniform grid.'
  end if

  return
end
subroutine nearest_point ( dim_num, c_num, center, point, nearest )
!
!*******************************************************************************
!
!! NEAREST_POINT finds the center point nearest a data point.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer C_NUM, the number of center points.
!
!    Input, real CENTER(DIM_NUM,C_NUM), the coordinates of the center points.
!
!    Input, real POINT(DIM_NUM), the data point to be checked.
!
!    Output, integer NEAREST, the index of the center point closest to
!    the data point.
!
  implicit none
!
  integer c_num
  integer dim_num
!
  real dist
  real dist_new
  real, dimension ( dim_num, c_num ) :: center
  integer j
  integer nearest
  real, dimension ( dim_num ) :: point
!
  dist = huge ( dist )
  nearest = 0

  do j = 1, c_num

    dist_new = sum ( ( point(1:dim_num) - center(1:dim_num,j) )**2 )

    if ( dist_new < dist ) then
      dist = dist_new
      nearest = j
    end if

  end do

  return
end
subroutine point_generate ( point_dist, dim_num, r_min, r_max, point_num, &
  point )
!
!*******************************************************************************
!
!! POINT_GENERATE generates data points for the problem.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_DIST, the point distribution to use.
!    1, use a uniform random distribution.
!    2, use a uniform grid of points.  (This hasn't been set up properly
!       except for 1D!).
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, real R_MIN(DIM_NUM), R_MAX(DIM_NUM), the coordinates of the
!    minimum and maximum corners of the region.
!
!    Input, integer POINT_NUM, the number of points to generate.
!
!    Output, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
  implicit none
!
  integer dim_num
  integer point_num
!
  integer i
  real point(dim_num,point_num)
  integer point_dist
  real r
  real r_max(dim_num)
  real r_min(dim_num)
!
  if ( point_dist == 1 ) then

    call random_number ( harvest = point(1:dim_num,1:point_num) )

    do i = 1, dim_num
      point(i,1:point_num) = r_min(i) + &
        point(i,1:point_num) * ( r_max(i) - r_min(i) )
    end do

  else if ( point_dist == 2 ) then

    do i = 1, point_num

      if ( point_num > 1 ) then
        r = real ( i - 1 ) / ( point_num - 1 )
      else
        r = 0.5E+00
      end if

      point(1:dim_num,i) = ( 1.0E+00 - r ) * r_min(1:dim_num) &
                                   + r   * r_max(1:dim_num)
    end do

  else

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POINT_GENERATE - Fatal error!'
    write ( *, '(a)' ) '  Meaningless input value of point distribution,'
    write ( *, '(a,i6)' ) '  POINT_DIST = ', point_dist
    stop

  end if

  return
end
subroutine point_print ( dim_num, point_num, point )
!
!*******************************************************************************
!
!! POINT_PRINT prints out the values of the data points.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
  implicit none
!
  integer dim_num
  integer point_num
!
  integer i
  real point(dim_num,point_num)
!
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Data points:'
  write ( *, '(a)' ) ' '

  do i = 1, point_num
    write ( *, '(i6,6f10.4)' ) i, point(1:dim_num,i)
  end do

  return
end
subroutine r_random ( rlo, rhi, r )
!
!*******************************************************************************
!
!! R_RANDOM returns a random real in a given range.
!
!
!  Modified:
!
!    06 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real RLO, RHI, the minimum and maximum values.
!
!    Output, real R, the randomly chosen value.
!
  implicit none
!
  real r
  real rhi
  real rlo
  real t
!
!  Pick T, a random number in (0,1).
!
  call random_number ( harvest = t )
!
!  Set R in ( RLO, RHI ).
!
  r = ( 1.0E+00 - t ) * rlo + t * rhi

  return
end
subroutine random_initialize ( seed )
!
!*******************************************************************************
!
!! RANDOM_INITIALIZE initializes the FORTRAN 90 random number seed.
!
!
!  Discussion:
!
!    If you don't initialize the random number generator, its behavior
!    is not specified.  If you initialize it simply by:
!
!      call random_seed
!
!    its behavior is not specified.  On the DEC ALPHA, if that's all you
!    do, the same random number sequence is returned.  In order to actually
!    try to scramble up the random number generator a bit, this routine
!    goes through the tedious process of getting the size of the random
!    number seed, making up values based on the current time, and setting
!    the random number seed.
!
!  Modified:
!
!    03 April 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, integer SEED.
!    If SEED is zero on input, then you're asking this routine to come up
!    with a seed value, which is returned as output.
!    If SEED is nonzero on input, then you're asking this routine to
!    use the input value of SEED to initialize the random number generator.
!
  implicit none
!
  integer count
  integer count_max
  integer count_rate
  integer i
  integer seed
  integer, allocatable :: seed_vector(:)
  integer seed_size
  real t
!
!  Initialize the random number seed.
!
  call random_seed
!
!  Determine the size of the random number seed.
!
  call random_seed ( size = seed_size )
!
!  Allocate a seed of the right size.
!
  allocate ( seed_vector(seed_size) )

  if ( seed /= 0 ) then

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RANDOM_INITIALIZE'
    write ( *, '(a,i12)' ) '  Initialize RANDOM_NUMBER with user SEED = ', seed

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'RANDOM_INITIALIZE'
    write ( *, '(a,i12)' ) '  Initialize RANDOM_NUMBER with arbitrary SEED = ', seed

  end if
!
!  Now set the seed.
!
  seed_vector(1:seed_size) = seed

  call random_seed ( put = seed_vector(1:seed_size) )
!
!  Free up the seed space.
!
  deallocate ( seed_vector )
!
!  Call the random number routine a bunch of times.
!
  do i = 1, 100
    call random_number ( harvest = t )
  end do

  return
end
subroutine rvec_range_input ( string, dim_num, value1, value2, ierror )
!
!*******************************************************************************
!
!! RVEC_RANGE_INPUT reads two real vectors from the user, representing a range.
!
!
!  Modified:
!
!    09 August 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) STRING, the prompt string.
!
!    Input, integer DIM_NUM, the number of dimensions.
!
!    Output, real VALUE1(DIM_NUM), VALUE2(DIM_NUM), the values entered 
!    by the user.
!
!    Output, integer IERROR, an error flag, which is zero if no error occurred.
!
  implicit none
!
  integer dim_num
!
  integer ierror
  character ( len = * ) string
  real value1(dim_num)
  real value2(dim_num)
!
  ierror = 0

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( string )
  read ( *, *, iostat = ierror ) value1(1:dim_num), value2(1:dim_num)

  return
end
subroutine timestamp ( )
!
!*******************************************************************************
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none
!
  character ( len = 8 ) ampm
  integer d
  character ( len = 8 ) date
  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
  character ( len = 10 )  time
  integer values(8)
  integer y
  character ( len = 5 ) zone
!
  call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end