subroutine ch_cap ( c )
!
!*******************************************************************************
!
!! CH_CAP capitalizes a single character.
!
!
!  Modified:
!
!    19 July 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character C, the character to capitalize.
!
  implicit none
!
  character c
  integer itemp
!
  itemp = ichar ( c )

  if ( 97 <= itemp .and. itemp <= 122 ) then
    c = char ( itemp - 32 )
  end if

  return
end
function ch_eqi ( c1, c2 )
!
!*******************************************************************************
!
!! CH_EQI is a case insensitive comparison of two characters for equality.
!
!
!  Examples:
!
!    CH_EQI ( 'A', 'a' ) is .TRUE.
!
!  Modified:
!
!    28 July 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C1, C2, the characters to compare.
!
!    Output, logical CH_EQI, the result of the comparison.
!
  implicit none
!
  logical ch_eqi
  character c1
  character c1_cap
  character c2
  character c2_cap
!
  c1_cap = c1
  c2_cap = c2

  call ch_cap ( c1_cap )
  call ch_cap ( c2_cap )

  if ( c1_cap == c2_cap ) then
    ch_eqi = .true.
  else
    ch_eqi = .false.
  end if

  return
end
subroutine ch_to_digit ( c, digit )
!
!*******************************************************************************
!
!! CH_TO_DIGIT returns the integer value of a base 10 digit.
!
!
!  Example:
!
!     C   DIGIT
!    ---  -----
!    '0'    0
!    '1'    1
!    ...  ...
!    '9'    9
!    ' '    0
!    'X'   -1
!
!  Modified:
!
!    04 August 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character C, the decimal digit, '0' through '9' or blank
!    are legal.
!
!    Output, integer DIGIT, the corresponding integer value.  If C was
!    'illegal', then DIGIT is -1.
!
  implicit none
!
  character c
  integer digit
!
  if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

    digit = ichar ( c ) - 48

  else if ( c == ' ' ) then

    digit = 0

  else

    digit = -1

  end if

  return
end
subroutine cluster_energy_compute ( dim_num, point_num, cluster_num, point, &
  cluster, cluster_center, cluster_energy )
!
!*******************************************************************************
!
!! CLUSTER_ENERGY_COMPUTE computes the energy of the clusters.
!
!
!  Modified:
!
!    18 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Output, integer CLUSTER(POINT_NUM), the cluster to which each
!    data point belongs.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the 
!    centers associated with the minimal energy clustering.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the energy associated with 
!    each cluster.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer, dimension ( point_num ) :: cluster
  real, dimension ( dim_num, cluster_num ) :: cluster_center
  real, dimension ( cluster_num ) :: cluster_energy
  integer i
  integer j
  real, dimension ( dim_num, point_num ) :: point
  real point_energy
!
  cluster_energy(1:cluster_num) = 0.0E+00

  do i = 1, point_num

    j = cluster(i)

    point_energy = sum ( &
      ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )

    cluster_energy(j) = cluster_energy(j) + point_energy

  end do

  return
end
subroutine cluster_initialize_1 ( dim_num, point_num, cluster_num, point, &
  cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_INITIALIZE_1 initializes the clusters to data points.
!
!
!  Discussion:
!
!    The cluster centers are simply chosen to be the first data points.
!
!  Modified:
!
!    06 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real cluster_center(dim_num,cluster_num)
  integer i
  real point(dim_num,point_num)
!
  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = point(1:dim_num,i)
  end do

  return
end
subroutine cluster_initialize_2 ( dim_num, point_num, cluster_num, point, &
  cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_INITIALIZE_2 initializes the cluster centers to random values.
!
!
!  Discussion:
!
!    In this case, the hyperbox containing the data is computed.
!
!    Then the cluster centers are chosen uniformly at random within
!    this hyperbox.
!
!    Of course, if the data is not smoothly distributed throughout
!    the box, many cluster centers will be isolated.
!
!  Modified:
!
!    05 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real cluster_center(dim_num,cluster_num)
  integer i
  real point(dim_num,point_num)
  real r(dim_num)
  real r_max(dim_num)
  real r_min(dim_num)
!
  r_min = minval ( point, 2 )
  r_max = maxval ( point, 2 )

  do i = 1, cluster_num

    call random_number ( harvest = r(1:dim_num) )

    cluster_center(1:dim_num,i) = &
      ( 1.0E+00 - r(1:dim_num) ) * r_min(1:dim_num) &
                + r(1:dim_num)   * r_max(1:dim_num)
  end do
    
  return
end
subroutine cluster_initialize_3 ( dim_num, point_num, cluster_num, point, &
  cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_INITIALIZE_3 initializes the cluster centers to random values.
!
!
!  Discussion:
!
!    In this case, each point is randomly assigned to a cluster, and
!    the cluster centers are then computed as the centroids of the points 
!    in the cluster.
!
!  Modified:
!
!    06 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real cluster_center(dim_num,cluster_num)
  integer cluster_population(cluster_num)
  integer i
  integer j
  real point(dim_num,point_num)
!
!  Assign one point to each cluster center.
!
  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = point(1:dim_num,i)
  end do

  cluster_population(1:cluster_num) = 1
!
!  The rest of the points get assigned randomly.
!
  do i = cluster_num+1, point_num
    call i_random ( 1, cluster_num, j )
    cluster_center(1:dim_num,j) = cluster_center(1:dim_num,j) + &
      point(1:dim_num,i)
    cluster_population(j) = cluster_population(j) + 1
  end do
!
!  Now average the points to get the centroid.
!
  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) / &
      real ( cluster_population(i) )
  end do
    
  return
end
subroutine cluster_initialize_4 ( dim_num, point_num, cluster_num, point, &
  cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_INITIALIZE_4 initializes the cluster centers to random values.
!
!
!  Discussion:
!
!    In this case, each data point is divided randomly among the
!    the cluster centers.
!
!  Modified:
!
!    06 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real cluster_center(dim_num,cluster_num)
  real cluster_factor(cluster_num)
  real cluster_weight(cluster_num)
  real divisor
  integer i
  integer j
  real point(dim_num,point_num)
!
  cluster_center(1:dim_num,1:cluster_num) = 0.0E+00
  cluster_weight(1:cluster_num) = 0.0E+00

  do i = 1, point_num

    call random_number ( harvest = cluster_factor(1:cluster_num) )
    divisor = sum ( cluster_factor(1:cluster_num) )
    cluster_factor(1:cluster_num) = cluster_factor(1:cluster_num) / divisor

    do j = 1, cluster_num
      cluster_center(1:dim_num,j) = cluster_center(1:dim_num,j) &
        + cluster_factor(j) * point(1:dim_num,i)
    end do

    cluster_weight(1:cluster_num) = cluster_weight(1:cluster_num) &
      + cluster_factor(1:cluster_num)

  end do
!
!  Now normalize,  so that each cluster center is now a convex 
!  combination of the points.
!
  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) &
      / cluster_weight(i)
  end do

  return
end
subroutine cluster_initialize_5 ( dim_num, point_num, cluster_num, point, &
  cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_INITIALIZE_5 initializes the cluster centers to random values.
!
!
!  Discussion:
!
!    In this case, each cluster center is a random convex combination 
!    of the data points.
!
!  Modified:
!
!    06 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real cluster_center(dim_num,cluster_num)
  real column_sum
  real factor(point_num,cluster_num)
  integer j
  real point(dim_num,point_num)
!
!  Get a PxC block of random factors.
!
  call random_number ( harvest = factor(1:point_num,1:cluster_num) )
!
!  Make each column of factors have unit sum.
!
  do j = 1, cluster_num
    column_sum = sum ( factor(1:point_num,j) )
    factor(1:point_num,j) = factor(1:point_num,j) / column_sum
  end do
!
!  Set centers = points * factors.
!
  cluster_center(1:dim_num,1:cluster_num) = &
    matmul ( point(1:dim_num,1:point_num), factor(1:point_num,1:cluster_num) )

  return
end
subroutine cluster_print_summary ( point_num, cluster_num, &
  cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! CLUSTER_PRINT_SUMMARY prints a summary of data about a clustering.
!
!
!  Modified:
!
!    05 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of
!    points assigned to each cluster.
!
!    Input, real CLUSTER_ENERGY(CLUSTER_NUM), the energy of the clusters.
!
  integer cluster_num
!
  real ce
  integer cep
  real ce_total
  real cluster_energy(cluster_num)
  integer cluster_population(cluster_num)
  integer cp
  integer cpp
  integer i
  integer point_num
!
  ce_total = sum ( cluster_energy(1:cluster_num) )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Clustering statistics:'
  write ( *, '(a)' ) ' '
  write ( *, '(a,i6)' ) '  Number of clusters is ', cluster_num
  write ( *, '(a,i6)' ) '  Number of points is   ', point_num
  write ( *, '(a,g14.6)' ) '  Total energy is       ', ce_total
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'Cluster   Population        Energy'
  write ( *, '(a)' ) '-------  -----------  -----------------'
  write ( *, '(a)' ) '              #    %      value       %'
  write ( *, '(a)' ) ' '

  do i = 1, cluster_num
    cp = cluster_population(i)
    cpp = int ( real ( 100 * cp ) / real ( point_num ) )
    ce = cluster_energy(i)
    cep = int ( ( ce * 100.0E+00 ) / ce_total )
    write ( *, '(4x,i3,2x,i6,2x,i3,g14.6,2x,i3)' ) i, cp, cpp, ce, cep
  end do

  return
end
subroutine cluster_write ( file_out_name, dim_num, point_num, cluster_num, &
  point, cluster, cluster_center )
!
!*******************************************************************************
!
!! CLUSTER_WRITE writes points and cluster centers to a file.
!
!
!  Modified:
!
!    09 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_OUT_NAME, the name of the output file.
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the points.
!
!    Input, integer CLUSTER(POINT_NUM), the cluster to which each point
!    belongs.
!
!    Input, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the cluster centers.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer cluster(point_num)
  real cluster_center(dim_num,cluster_num)
  character ( len = * ) file_out_name
  integer file_out_unit
  integer i
  integer ierror
  integer ios
  real point(dim_num,point_num)
!
  call get_unit ( file_out_unit )

  open ( unit = file_out_unit, file = file_out_name, status = 'replace', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'CLUSTER_WRITE - Fatal error!'
    write ( *, '(a)' ) '  Could not open the output file: ' // &
      trim ( file_out_name )
    stop
  end if

  write ( file_out_unit, '(a)' ) '# ' // trim ( file_out_name )
  write ( file_out_unit, '(a)' ) '#  Created by CLUSTER_WRITE'
  write ( file_out_unit, '(a)' ) '#'
  write ( file_out_unit, '(a)' ) '#  Results of a K-Means computation.'
  write ( file_out_unit, '(a)' ) '#'
  write ( file_out_unit, '(a)' ) '#  First we list the points.'
  write ( file_out_unit, '(a)' ) '#  Then the cluster centers.'
  write ( file_out_unit, '(a)' ) '#  Cluster centers are set off by blanks.'
  write ( file_out_unit, '(a)' ) '#'
  write ( file_out_unit, '(a)' ) '#  Points and centers list their cluster.'
  write ( file_out_unit, '(a)' ) '#'

  do i = 1, point_num
    write ( file_out_unit, '(2g14.6,2x,i3)' ) point(1:dim_num,i), cluster(i)
  end do

  do i = 1, cluster_num
    write ( file_out_unit, '(a)' ) ' '
    write ( file_out_unit, '(2g14.6,2x,i3)' ) cluster_center(1:dim_num,i), i
  end do

  close ( unit = file_out_unit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'CLUSTER_WRITE:'
  write ( *, '(a,i6)' ) '  Wrote cluster centers to file.'

  return
end
subroutine file_column_count ( file_name, ncolumn )
!
!*******************************************************************************
!
!! FILE_COLUMN_COUNT counts the number of columns in the first line of a file.
!
!
!  Discussion:
!
!    The file is assumed to be a simple text file.
!
!    Most lines of the file is presumed to consist of NCOLUMN words, separated
!    by spaces.  There may also be some blank lines, and some comment lines,
!    which have a "#" in column 1.
!
!    The routine tries to find the first non-comment non-blank line and
!    counts the number of words in that line.
!
!    If all lines are blanks or comments, it goes back and tries to analyze
!    a comment line.
!
!  Modified:
!
!    21 June 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_NAME, the name of the file.
!
!    Output, integer NCOLUMN, the number of columns assumed to be in the file.
!
  implicit none
!
  character ( len = * ) file_name
  logical got_one
  integer ios
  integer iunit
  character ( len = 256 ) line
  integer ncolumn
!
!  Open the file.
!
  call get_unit ( iunit )

  open ( unit = iunit, file = file_name, status = 'old', form = 'formatted', &
    access = 'sequential', iostat = ios )

  if ( ios /= 0 ) then
    ncolumn = - 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the file:'
    write ( *, '(a)' ) '    ' // trim ( file_name )
    return
  end if
!
!  Read one line, but skip blank lines and comment lines.
!
  got_one = .false.

  do

    read ( iunit, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      exit
    end if

    if ( len_trim ( line ) == 0 ) then
      cycle
    end if

    if ( line(1:1) == '#' ) then
      cycle
    end if

    got_one = .true.
    exit

  end do

  if ( .not. got_one ) then

    rewind ( iunit )

    do

      read ( iunit, '(a)', iostat = ios ) line

      if ( ios /= 0 ) then
        exit
      end if

      if ( len_trim ( line ) == 0 ) then
        cycle
      end if

      got_one = .true.
      exit

    end do

  end if

  close ( unit = iunit )

  if ( .not. got_one ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!'
    write ( *, '(a)' ) '  The file does not seem to contain any data.'
    ncolumn = 0
    return
  end if

  call word_count ( line, ncolumn )

  return
end
subroutine get_unit ( iunit )
!
!*******************************************************************************
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!  Modified:
!
!    02 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5 and 6).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
  implicit none
!
  integer i
  integer ios
  integer iunit
  logical lopen

  iunit = 0

  do i = 1, 99

    if ( i /= 5 .and. i /= 6 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )

      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if

  end do

  return
end
subroutine hmeans_01 ( dim_num, point_num, cluster_num, it_max, it, point, &
  cluster, cluster_center, cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! HMEANS_01 applies the H-Means algorithm to a partition problem.
!
!
!  Discussion:
!
!    The data for the H-Means problem is a set of N points X in
!    M-dimensions, and a desired number of clusters K.
!
!    The goal is to determine K points Z, called cluster centers, so that
!    if we associate each point X with its nearest Z value, we minimize
!    the standard deviation or cluster energy.  Writing CLUSTER(I) to
!    indicate the index of the nearest cluster center to point X(I), the
!    energy can be written as:
!
!      Energy = Sum ( 1 <= I <= N ) || X(I) - Z(CLUSTER(I)) ||**2
!
!    where
!
!      || X - Z ||**2 = Sum ( 1 <= J <= M ) ( X(J) - Z(J) )**2
!
!  Reference:
!
!    Wendy Martinez and Angel Martinez,
!    Computational Statistics Handbook with MATLAB,
!    pages 373-376,
!    Chapman and Hall / CRC, 2002.
!
!  Modified:
!
!    22 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer IT_MAX, the maximum number of iterations.
!
!    Output, integer IT, the number of iterations taken.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Output, integer CLUSTER(POINT_NUM), the cluster to which each
!    data point belongs.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the 
!    centers associated with the minimal energy clustering.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the energy associated with 
!    each cluster.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer c
  real, dimension ( dim_num, cluster_num ) :: centroid
  integer, dimension ( point_num ) :: cluster
  real, dimension ( dim_num, cluster_num ) :: cluster_center
  real, dimension ( cluster_num ) :: cluster_energy
  integer, dimension ( cluster_num ) :: cluster_population
  logical, parameter :: debug = .true.
  integer i
  integer it
  integer it_max
  integer j
  integer missed
  real, dimension ( dim_num, point_num ) :: point
  real point_energy
  real point_energy_min
  integer swap
!
  do it = 1, it_max
!
!  #1:
!  Assign each point to the cluster of its nearest center.
!
    swap = 0

    do i = 1, point_num

      point_energy_min = huge ( point_energy_min )
      c = cluster(i)

      do j = 1, cluster_num

        point_energy = sum ( &
          ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )

        if ( point_energy < point_energy_min ) then
          point_energy_min = point_energy
          cluster(i) = j
        end if

      end do

      if ( c /= cluster(i) ) then
        swap = swap + 1
      end if
      
    end do
!
!  Terminate if no points were swapped.
!
    if ( it > 1 ) then
      if ( swap == 0 ) then
        exit
      end if
    end if
!
!  #2:
!  Determine the total energy of the new clustering with current centroids.
!
    call cluster_energy_compute ( dim_num, point_num, cluster_num, point, &
      cluster, cluster_center, cluster_energy )
!
!  #3:
!  Determine the centroids of the clusters.
!
    centroid(1:dim_num,1:cluster_num) = 0.0E+00
    cluster_population(1:cluster_num) = 0

    do i = 1, point_num
      j = cluster(i)
      cluster_population(j) = cluster_population(j) + 1
      centroid(1:dim_num,j) = centroid(1:dim_num,j) + point(1:dim_num,i)
    end do
!
!  Now divide by the population to get the centroid.
!  But if a center has no population, pick a point at random.
!
    missed = 0

    do j = 1, cluster_num

      if ( cluster_population(j) /= 0 ) then
        centroid(1:dim_num,j) = centroid(1:dim_num,j) &
          / real ( cluster_population(j) )
      else
        missed = missed + 1

        centroid(1:dim_num,j) = point(1:dim_num,missed)

      end if

    end do

    cluster_center(1:dim_num,1:cluster_num) = centroid(1:dim_num,1:cluster_num)
!
!  #4:
!  Determine the total energy of the current clustering with new centroids.
!
    call cluster_energy_compute ( dim_num, point_num, cluster_num, point, &
      cluster, cluster_center, cluster_energy )
 
  end do

  return
end
subroutine hmeans_02 ( dim_num, point_num, cluster_num, it_max, it, point, &
  cluster, cluster_center, cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! HMEANS_02 applies the H-Means algorithm to a partition problem.
!
!
!  Discussion:
!
!    This is a simple routine to group a set of points into K clusters,
!    each with a center point, in such a way that the total cluster 
!    energy is minimized.  The total cluster energy is the sum of the
!    squares of the distances of each point to the center of its cluster.
!
!    The algorithm begins with an initial estimate for the cluster centers:
!
!    1. The points are assigned to the nearest cluster centers.
!
!    2. The iteration stops if the total energy has not changed 
!        significantly, or we have reached the maximum number of iterations.
!
!    3. Each cluster center is replaced by the centroid of the points
!       in the cluster.
!
!    4. Return to step 1.
!
!    The algorithm may fail to find the best solution.
!
!  Modified:
!
!    22 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer IT_MAX, the maximum number of iterations to carry out.
!
!    Output, integer IT, the number of iterations taken.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, integer CLUSTER(POINT_NUM), the cluster to which each
!    point belongs.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the coordinates
!    of the cluster centers.
!
!    Output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of
!    points assigned to each cluster.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the energy of the clusters.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer c
  integer cluster(point_num)
  real cluster_center(dim_num,cluster_num)
  real cluster_energy(cluster_num)
  integer cluster_population(cluster_num)
  logical, parameter :: debug = .false.
  real energy(cluster_num)
  integer i
  integer it
  integer it_max
  integer j
  integer list(1)
  real point(dim_num,point_num)
  integer swap
!
!  Idiot checks.
!  (You laugh, but wait til you've spent an hour tracking down
!  bugs caused by a variable that was never set.)
!
  if ( cluster_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  CLUSTER_NUM < 1.'
    stop
  end if

  if ( dim_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  DIM_NUM < 1.'
    stop
  end if

  if ( point_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  POINT_NUM < 1.'
    stop
  end if

  if ( it_max < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'HMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  IT_MAX < 0.'
    stop
  end if

  it = 0

  do
!
!  Given centers, assign points to nearest center.
!
    cluster_population(1:cluster_num) = 0
    cluster_energy(1:cluster_num) = 0.0E+00

    swap = 0

    do i = 1, point_num

      do j = 1, cluster_num
        energy(j) = sum ( &
          ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )
      end do

      list = minloc ( energy(1:cluster_num) )

      c = list(1)

      if ( c /= cluster(i) ) then
        swap = swap + 1
      end if

      cluster(i) = c
      cluster_energy(c) = cluster_energy(c) + energy(c) 
      cluster_population(c) = cluster_population(c) + 1

    end do

    if ( debug ) then
      write ( *, '(i3,g14.6)' ) it, sum ( cluster_energy(1:cluster_num) )
    end if

    if ( it > 0 ) then
      if ( swap == 0 ) then
        exit
      end if
    end if

    if ( it >= it_max ) then
      exit
    end if

    it = it + 1
!
!  Given points in cluster, replace center by centroid.
!
    cluster_center(1:dim_num,1:cluster_num) = 0.0E+00
  
    do i = 1, point_num
      c = cluster(i)
      cluster_center(1:dim_num,c) = cluster_center(1:dim_num,c) &
        + point(1:dim_num,i)
    end do

    do i = 1, cluster_num
      if ( cluster_population(i) /= 0 ) then
        cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) / &
          real ( cluster_population(i) )
      else
        call i_random ( 1, point_num, j )
        cluster_center(1:dim_num,i) = point(1:dim_num,j)
      end if
    end do

  end do

  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 kmeans_01 ( dim_num, point_num, cluster_num, it_max, it, point, &
  cluster, cluster_center, cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! KMEANS_01 applies the K-Means algorithm to a partition problem.
!
!
!  Discussion:
!
!    Given a matrix of POINT_NUM observations on DIM_NUM variables, the
!    observations are to be allocated to CLUSTER_NUM clusters in such 
!    a way that the within-cluster sum of squares is minimized.
!
!  Reference:
!
!    D N Sparks,
!    Algorithm AS 58: Euclidean Cluster Analysis,
!    Applied Statistics,
!    Volume 22, Number 1, 1973, pages 126-130.
!
!  Modified:
!
!    23 March 2002
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer IT_MAX, the maximum number of iterations to carry out.
!
!    Output, integer IT, the number of iterations taken.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the points.
!
!    Output, integer CLUSTER(POINT_NUM), indicates which cluster
!    each point belongs to.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the cluster
!    centers.
!
!    Output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of points
!    in each cluster.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the cluster energies.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer cluster(point_num)
  real cluster_center(dim_num,cluster_num)
  real cluster_energy(cluster_num)
  integer cluster_population(cluster_num)
  real dc
  real de
  real f(point_num)
  integer i
  integer il
  integer ir
  integer it
  integer it_max
  integer j
  integer k
  integer list(1)
  real point(dim_num,point_num)
  integer swap
!
  it = 0
!
!  Idiot checks.
!
  if ( cluster_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_01 - Fatal error!'
    write ( *, '(a)' ) '  CLUSTER_NUM < 1.'
    stop
  end if

  if ( dim_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_01 - Fatal error!'
    write ( *, '(a)' ) '  DIM_NUM < 1.'
    stop
  end if

  if ( point_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_01 - Fatal error!'
    write ( *, '(a)' ) '  POINT_NUM < 1.'
    stop
  end if
!
!  For each observation, calculate the distance from each cluster
!  center, and assign to the nearest.
!
  do i = 1, point_num

    do j = 1, cluster_num
      cluster_energy(j) = sum ( &
        ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )
    end do

    list = minloc ( cluster_energy(1:cluster_num) )
    cluster(i) = list(1)

  end do
!
!  Determine the cluster population counts.
!
  cluster_population(1:cluster_num) = 0

  do i = 1, point_num
    j = cluster(i)
    cluster_population(j) = cluster_population(j) + 1
  end do
!
!  Calculate the mean and sum of squares for each cluster.
!
  cluster_center(1:dim_num,1:cluster_num) = 0.0E+00

  do i = 1, point_num
    j = cluster(i)
    cluster_center(1:dim_num,j) = cluster_center(1:dim_num,j) &
      + point(1:dim_num,i)
  end do

  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) &
      / real ( cluster_population(i) )
  end do
!
!  Set the point energies.
!
  f(1:point_num) = 0.0E+00

  do i = 1, point_num
    j = cluster(i)
    f(i) = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )
  end do
!
!  Set the cluster energies.
!
  cluster_energy(1:cluster_num) = 0.0E+00

  do i = 1, point_num
    j = cluster(i)
    cluster_energy(j) = cluster_energy(j) + f(i)
  end do
!
!  Adjust the point energies by a weight factor.
!
  do i = 1, point_num
    j = cluster(i)
    if ( cluster_population(j) > 1 ) then
      f(i) = f(i) * real ( cluster_population(j) ) &
        / real ( cluster_population(j) - 1 )
    end if
  end do
!
!  Examine each observation in turn to see if it should be
!  reassigned to a different cluster.
!
  do

    if ( it >= it_max ) then
      exit
    end if

    it = it + 1

    swap = 0

    do i = 1, point_num

      il = cluster(i)
      ir = il

      if ( cluster_population(il) <= 1 ) then
        cycle
      end if

      dc = f(i)

      do j = 1, cluster_num

        if ( j /= il ) then

          de = sum ( &
            ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 ) &
            * real ( cluster_population(j) ) &
            / real ( cluster_population(j) + 1 )

          if ( de < dc ) then
            dc = de
            ir = j
          end if

        end if

      end do
!
!  If the lowest value was obtained by staying in the current cluster,
!  then cycle.
!
      if ( ir == il ) then
        cycle
      end if
!
!  Reassign the point from cluster IL to cluster IR.
!
      cluster_center(1:dim_num,il) = &
        ( cluster_center(1:dim_num,il) * real ( cluster_population(il) ) &
        - point(1:dim_num,i) ) / real ( cluster_population(il) - 1 )

      cluster_center(1:dim_num,ir) = &
        ( cluster_center(1:dim_num,ir) * real ( cluster_population(ir) ) &
        + point(1:dim_num,i) ) / real ( cluster_population(ir) + 1 )

      cluster_energy(il) = cluster_energy(il) - f(i)
      cluster_energy(ir) = cluster_energy(ir) + dc
      cluster_population(ir) = cluster_population(ir) + 1
      cluster_population(il) = cluster_population(il) - 1

      cluster(i) = ir
!
!  Adjust the value of F for points in clusters IL and IR.
!
      do j = 1, point_num

        k = cluster(j)

        if ( k == il .or. k == ir ) then

          f(j) = sum ( &
           ( point(1:dim_num,j) - cluster_center(1:dim_num,k) )**2 )

          f(j) = f(j) * real ( cluster_population(k) ) &
            / ( real ( cluster_population(k) - 1 ) )

        end if

      end do

      swap = swap + 1

    end do
!
!  Exit if no reassignments were made during this iteration.
!
    if ( swap == 0 ) then
      exit
    end if

  end do

  return
end
subroutine kmeans_02 ( dim_num, point_num, cluster_num, it_max, it, point, &
  cluster, cluster_center, cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! KMEANS_02 applies the K-Means algorithm to a partition problem.
!
!
!  Discussion:
!
!    The routine attempts to divide POINT_NUM points in 
!    DIM_NUM-dimensional space into CLUSTER_NUM clusters so that the within 
!    cluster sum of squares is minimized.
!
!  Reference:
!
!    J A Hartigan, M A Wong,
!    Algorithm AS 136: A K-Means Clustering Algorithm,
!    Applied Statistics,
!    Volume 28, Number 1, 1979, pages 100-108.
!
!  Modified:
!
!    10 March 2002
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer IT_MAX, the maximum number of iterations allowed.
!
!    Output, integer IT, the number of iterations taken.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Output, integer CLUSTER(POINT_NUM), the cluster each point belongs to.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the cluster
!    centers.
!
!    Output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of points 
!    in each cluster.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the within-cluster sum 
!    of squares.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real an1(cluster_num)
  real an2(cluster_num)
  integer cluster(point_num)
  integer cluster2(point_num)
  real cluster_center(dim_num,cluster_num)
  real cluster_energy(cluster_num)
  integer cluster_population(cluster_num)
  real d(point_num)
  real dc
  real db
  real dt(2)
  integer i 
  integer ifault
  integer il
  integer indx
  integer it
  integer it2
  integer it_max
  integer itran(cluster_num)
  integer j
  integer l
  integer live(cluster_num)
  integer ncp(cluster_num)
  real point(dim_num,point_num)
  real temp
!
  it = 0

  if ( cluster_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  CLUSTER_NUM < 1.'
    stop
  end if

  if ( cluster_num >= point_num ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  CLUSTER_NUM >= POINT_NUM.'
    stop
  end if

  if ( dim_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  DIM_NUM < 1.'
    stop
  end if

  if ( point_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_02 - Fatal error!'
    write ( *, '(a)' ) '  POINT_NUM < 1.'
    stop
  end if
!
!  For each point I, find its two closest centers, CLUSTER(I) and CLUSTER2(I).
!  Assign it to CLUSTER(I).
!
  do i = 1, point_num

    cluster(i) = 1
    cluster2(i) = 2

    do il = 1, 2
      dt(il) = sum ( &
        ( point(1:dim_num,i) - cluster_center(1:dim_num,il) )**2 )
    end do

    if ( dt(1) > dt(2) ) then
      cluster(i) = 2
      cluster2(i) = 1
      temp = dt(1)
      dt(1) = dt(2)
      dt(2) = temp
    end if

    do l = 3, cluster_num

      db = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l) )**2 )

      if ( db < dt(1) ) then
        dt(2) = dt(1)
        cluster2(i) = cluster(i)
        dt(1) = db
        cluster(i) = l
      else if ( db < dt(2) ) then
        dt(2) = db
        cluster2(i) = l
      end if

    end do

  end do
!
!  Update cluster centers to be the average of points contained
!  within them.
!
  cluster_population(1:cluster_num) = 0
  cluster_center(1:dim_num,1:cluster_num) = 0.0E+00

  do i = 1, point_num
    l = cluster(i)
    cluster_population(l) = cluster_population(l) + 1
    cluster_center(1:dim_num,l) = cluster_center(1:dim_num,l) &
      + point(1:dim_num,i)
  end do
!
!  Check to see if there is any empty cluster.
!
  do l = 1, cluster_num

    if ( cluster_population(l) == 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'KMEANS_02 - Fatal error!'
      write ( *, '(a)' ) '  There is an empty cluster.'
      stop
    end if

    cluster_center(1:dim_num,l) = cluster_center(1:dim_num,l) &
      / real ( cluster_population(l) )
!
!  Initialize AN1, AN2, ITRAN and NCP
!  AN1(L) = CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) - 1)
!  AN2(L) = CLUSTER_POPULATION(L) / (CLUSTER_POPULATION(L) + 1)
!  ITRAN(L) = 1 if cluster L is updated in the quick-transfer stage,
!           = 0 otherwise
!  In the optimal-transfer stage, NCP(L) stores the step at which
!  cluster L is last updated.
!  In the quick-transfer stage, NCP(L) stores the step at which
!  cluster L is last updated plus POINT_NUM.
!
    an2(l) = real ( cluster_population(l) ) / real ( cluster_population(l) + 1 )

    if ( cluster_population(l) > 1 ) then
      an1(l) = real ( cluster_population(l) ) &
        / real ( cluster_population(l) - 1 )
    else
      an1(l) = huge ( an1(l) )
    end if

    itran(l) = 1
    ncp(l) = -1

  end do

  indx = 0
  ifault = 2

  do it2 = 1, it_max

    it = it2
!
!  In this stage, there is only one pass through the data.   Each
!  point is re-allocated, if necessary, to the cluster that will
!  induce the maximum reduction in within-cluster sum of squares.
!
    call kmeans_02_optra ( dim_num, point_num, cluster_num, point, &
      cluster_center, cluster, cluster2, cluster_population, an1, an2, &
      ncp, d, itran, live, indx )
!
!  Stop if no transfer took place in the last POINT_NUM optimal transfer steps.
!
    if ( indx == point_num ) then
      ifault = 0
      exit
    end if
!
!  Each point is tested in turn to see if it should be re-allocated
!  to the cluster to which it is most likely to be transferred,
!  CLUSTER2(I), from its present cluster, CLUSTER(I).   Loop through the
!  data until no further change is to take place.
!
    call kmeans_02_qtran ( dim_num, point_num, cluster_num, point, &
      cluster_center, cluster, cluster2, cluster_population, an1, an2, &
      ncp, d, itran, indx )
!
!  If there are only two clusters, there is no need to re-enter the
!  optimal transfer stage.
!
    if ( cluster_num == 2 ) then
      ifault = 0
      exit
    end if
!
!  NCP has to be set to 0 before entering OPTRA.
!
    ncp(1:cluster_num) = 0

  end do

  if ( ifault == 2 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_02 - Warning!'
    write ( *, '(a)' ) '  Maximum number of iterations reached'
    write ( *, '(a)' ) '  without convergence.'
  end if
!
!  Compute the within-cluster sum of squares for each cluster.
!
  cluster_center(1:dim_num,1:cluster_num) = 0.0E+00

  do i = 1, point_num
    cluster_center(1:dim_num,cluster(i)) = &
    cluster_center(1:dim_num,cluster(i)) + point(1:dim_num,i)
  end do

  do j = 1, dim_num
    cluster_center(j,1:cluster_num) = cluster_center(j,1:cluster_num) &
      / real ( cluster_population(1:cluster_num) )
  end do

  call cluster_energy_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_energy )

  return
end
subroutine kmeans_02_optra ( dim_num, point_num, cluster_num, point, &
  cluster_center, cluster, cluster2, cluster_population, an1, an2, ncp, &
  d, itran, live, indx )
!
!*******************************************************************************
!
!! KMEANS_02_OPTRA carries out the optimal transfer stage.
!
!
!  Discussion:
!
!    Each point is re-allocated, if necessary, to the cluster that
!    will induce a maximum reduction in the within-cluster sum of
!    squares.
!
!  Reference:
!
!    J A Hartigan, M A Wong,
!    A K-Means Clustering Algorithm,
!    Applied Statistics,
!    Volume 28, Number 1, 1979, pages 100-108.
!
!  Modified:
!
!    08 March 2002
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the cluster
!    centers.
!
!    Input/output, integer CLUSTER(POINT_NUM), the cluster each point 
!    belongs to.
!
!    Input/output, integer CLUSTER2(POINT_NUM), the cluster to which
!    each point is most likely to be transferred to.
!
!    Input/output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of 
!    points in each cluster.
!
!    ?, real AN1(CLUSTER_NUM), ?
!
!    ?, real AN2(CLUSTER_NUM), ?
!
!    ?, integer NCP(CLUSTER_NUM), ?
!
!    ?, real D(POINT_NUM), ?
!
!    Input/output, integer ITRAN(CLUSTER_NUM), ?
!
!    Input/output, integer LIVE(CLUSTER_NUM), ?
!
!    Input/output, integer INDX, ?
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real al1
  real al2
  real alt
  real alw
  real an1(cluster_num)
  real an2(cluster_num)
  integer cluster(point_num)
  integer cluster2(point_num)
  real cluster_center(dim_num,cluster_num)
  integer cluster_population(cluster_num)
  real d(point_num)
  real dc
  real dd
  integer i
  integer indx
  integer itran(cluster_num)
  integer j
  integer l
  integer l1
  integer l2
  integer live(cluster_num)
  integer ll
  integer ncp(cluster_num)
  real point(dim_num,point_num)
  real r2
  real rr
!
!  If cluster L is updated in the last quick-transfer stage, it
!  belongs to the live set throughout this stage.   Otherwise, at
!  each step, it is not in the live set if it has not been updated
!  in the last POINT_NUM optimal transfer steps.
!
  do l = 1, cluster_num
    if ( itran(l) == 1 ) then
      live(l) = point_num + 1
    end if
  end do

  do i = 1, point_num

    indx = indx + 1
    l1 = cluster(i)
    l2 = cluster2(i)
    ll = l2
!
!  If point I is the only member of cluster L1, no transfer.
!
    if ( cluster_population(l1) > 1 ) then
!
!  If L1 has been updated in this stage, re-compute D(I).
!
      if ( ncp(l1) /= 0 ) then
        d(i) = an1(l1) * sum ( &
          ( point(1:dim_num,i) - cluster_center(1:dim_num,l1) )**2 )
      end if
!
!  Find the cluster with minimum R2.
!
      r2 = an2(l2) * sum ( &
        ( point(1:dim_num,i) - cluster_center(1:dim_num,l2) )**2 )

      do l = 1, cluster_num
!
!  If I >= LIVE(L1), then L1 is not in the live set.   If this is
!  true, we only need to consider clusters that are in the live set
!  for possible transfer of point I.   
!
!  Otherwise, we need to consider all possible clusters.
!
!  BADLY FORMED LOGICAL HERE.
!
        if ( i >= live(l1) .and. i >= live(l) .or. l == l1 .or. l == ll ) then
          go to 60
        end if

        rr = r2 / an2(l)

        dc = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l) )**2 )

        if ( dc < rr ) then
          r2 = dc * an2(l)
          l2 = l
        end if

60      continue

      end do
!
!  If no transfer is necessary, L2 is the new CLUSTER2(I).
! 
      if ( r2 >= d(i) ) then

        cluster2(i) = l2

      else
!
!  Update cluster centers, LIVE, NCP, AN1 and AN2 for clusters L1 and
!  L2, and update CLUSTER(I) and CLUSTER2(I).
!
        indx = 0
        live(l1) = point_num + i
        live(l2) = point_num + i
        ncp(l1) = i
        ncp(l2) = i
        al1 = cluster_population(l1)
        alw = al1 - 1.0E+00
        al2 = cluster_population(l2)
        alt = al2 + 1.0E+00

        cluster_center(1:dim_num,l1) = ( cluster_center(1:dim_num,l1) * al1 &
          - point(1:dim_num,i) ) / alw

        cluster_center(1:dim_num,l2) = ( cluster_center(1:dim_num,l2) * al2 &
          + point(1:dim_num,i) ) / alt

        cluster_population(l1) = cluster_population(l1) - 1
        cluster_population(l2) = cluster_population(l2) + 1
        an2(l1) = alw / al1

        if ( alw > 1.0E+00 ) then
          an1(l1) = alw / ( alw - 1.0E+00 )
        else
          an1(l1) = huge ( an1(l1) )
        end if

        an1(l2) = alt / al2
        an2(l2) = alt / ( alt + 1.0E+00 )
        cluster(i) = l2
        cluster2(i) = l1

      end if

    end if

    if ( indx == point_num ) then
      return
    end if

  end do
!
!  ITRAN(L) = 0 before entering QTRAN.
!
  itran(1:cluster_num) = 0
!
!  LIVE(L) has to be decreased by POINT_NUM before re-entering OPTRA.
!
  live(1:cluster_num) = live(1:cluster_num) - point_num

  return
end
subroutine kmeans_02_qtran ( dim_num, point_num, cluster_num, point, &
  cluster_center, cluster, cluster2, cluster_population, an1, an2, ncp, &
  d, itran, indx )
!
!*******************************************************************************
!
!! KMEANS_02_QTRAN carries out the quick transfer stage.
!
!
!  Discussion:
!
!    For each point I, CLUSTER(I) and CLUSTER2(I) are switched, if necessary, 
!    to reduce within-cluster sum of squares.  The cluster centers are
!    updated after each step.
!
!  Reference:
!
!    J A Hartigan, M A Wong,
!    A K-Means Clustering Algorithm,
!    Applied Statistics,
!    Volume 28, Number 1, 1979, pages 100-108.
!
!  Modified:
!
!    10 March 2002
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the coordinates of the points.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the cluster
!    centers.
!
!    Input/output, integer CLUSTER(POINT_NUM), the cluster each point 
!    belongs to.
!
!    Input/output, integer CLUSTER2(POINT_NUM), the cluster to which
!    each point is most likely to be transferred to.
!
!    Input/output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number of 
!    points in each cluster.
!
!    ?, real AN1(CLUSTER_NUM), ?
!
!    ?, real AN2(CLUSTER_NUM), ?
!
!    ?, integer NCP(CLUSTER_NUM), ?
!
!    ?, real D(POINT_NUM), ?
!
!    ?, integer ITRAN(CLUSTER_NUM), ?
!
!    Input/output, integer INDX, is set to 0 if any updating occurs.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  real al1
  real al2
  real alt
  real alw
  real an1(cluster_num)
  real an2(cluster_num)
  integer cluster(point_num)
  integer cluster2(point_num)
  real cluster_center(dim_num,cluster_num)
  integer cluster_population(cluster_num)
  integer count
  real d(point_num)
  real dd
  integer i
  integer indx
  integer itran(cluster_num)
  integer j
  integer l1
  integer l2
  integer ncp(cluster_num)
  real point(dim_num,point_num)
  real r2
  integer step
!
!  In the optimal transfer stage, NCP(L) indicates the step at which
!  cluster L is last updated.   In the quick transfer stage, NCP(L)
!  is equal to the step at which cluster L is last updated plus POINT_NUM.
!
  count = 0
  step = 0

  do

    do i = 1, point_num

      count = count + 1
      step = step + 1
      l1 = cluster(i)
      l2 = cluster2(i)
!
!  If point I is the only member of cluster L1, no transfer.
!
      if ( cluster_population(l1) > 1 ) then
!
!  If STEP > NCP(L1), no need to re-compute distance from point I to
!  cluster L1.   Note that if cluster L1 is last updated exactly POINT_NUM
!  steps ago, we still need to compute the distance from point I to
!  cluster L1.
!
        if ( step <= ncp(l1) ) then
          d(i) = an1(l1) * sum ( &
            ( point(1:dim_num,i) - cluster_center(1:dim_num,l1) )**2 )
        end if
!
!  If STEP >= both NCP(L1) and NCP(L2) there will be no transfer of
!  point I at this step.
!
        if ( step < ncp(l1) .or. step < ncp(l2) ) then

          r2 = d(i) / an2(l2)

          dd = sum ( ( point(1:dim_num,i) - cluster_center(1:dim_num,l2) )**2 )
!
!  Update cluster centers, NCP, CLUSTER_POPULATION, ITRAN, AN1 and AN2 
!  for clusters L1 and L2.   Also update CLUSTER(I) and CLUSTER2(I).   
!
!  Note that if any updating occurs in this stage, INDX is set back to 0.
!
          if ( dd < r2 ) then

            count = 0
            indx = 0
            itran(l1) = 1
            itran(l2) = 1
            ncp(l1) = step + point_num
            ncp(l2) = step + point_num
            al1 = cluster_population(l1)
            alw = al1 - 1.0E+00
            al2 = cluster_population(l2)
            alt = al2 + 1.0E+00

            cluster_center(1:dim_num,l1) = &
              ( cluster_center(1:dim_num,l1) * al1 &
              - point(1:dim_num,i) ) / alw

            cluster_center(1:dim_num,l2) = &
              ( cluster_center(1:dim_num,l2) * al2 &
              + point(1:dim_num,i) ) / alt

            cluster_population(l1) = cluster_population(l1) - 1
            cluster_population(l2) = cluster_population(l2) + 1
            an2(l1) = alw / al1

            if ( alw > 1.0E+00 ) then
              an1(l1) = alw / ( alw - 1.0E+00 )
            else
              an1(l1) = huge ( an1(l1) )
            end if

            an1(l2) = alt / al2
            an2(l2) = alt / ( alt + 1.0E+00 )
            cluster(i) = l2
            cluster2(i) = l1

          end if

        end if

      end if
!
!  If no re-allocation took place in the last POINT_NUM steps, return.
!
      if ( count == point_num ) then
        return
      end if

    end do

  end do

  return
end
subroutine kmeans_03 ( dim_num, point_num, cluster_num, it_max, it, point, &
  cluster, cluster_center, cluster_population, cluster_energy )
!
!*******************************************************************************
!
!! KMEANS_03 applies the K-Means algorithm to a partition problem..
!
!
!  Discussion:
!
!    It is possible for a straightforward K-Means algorithm to
!    halt at a non-optimal partition of the points.  This routine
!    tries to improve the input partition if possible.
!
!  Reference:
!
!    Wendy Martinez and Angel Martinez,
!    Computational Statistics Handbook with MATLAB,
!    pages 373-376,
!    Chapman and Hall / CRC, 2002.
!
!  Modified:
!
!    22 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of data points.
!
!    Input, integer CLUSTER_NUM, the number of clusters.
!
!    Input, integer IT_MAX, the maximum number of iterations to carry out.
!
!    Output, integer IT, the number of iterations taken.
!
!    Input, real POINT(DIM_NUM,POINT_NUM), the data points.
!
!    Input/output, integer CLUSTER(POINT_NUM), the cluster to which
!    each point belongs.  On output, these may have been altered.
!
!    Input/output, real CLUSTER_CENTER(DIM_NUM,CLUSTER_NUM), the 
!    centers associated with the clustering.  On output, these may 
!    have been altered.
!
!    Output, integer CLUSTER_POPULATION(CLUSTER_NUM), the number
!    of points in each cluster.
!
!    Output, real CLUSTER_ENERGY(CLUSTER_NUM), the energy of the clusters.
!
  implicit none
!
  integer cluster_num
  integer dim_num
  integer point_num
!
  integer ci
  integer cj
  integer, dimension ( point_num ) :: cluster
  real, dimension ( dim_num, cluster_num ) :: cluster_center
  real, dimension ( cluster_num ) :: cluster_energy
  integer, dimension ( cluster_num ) :: cluster_population
  logical, parameter :: debug = .true.
  real, dimension ( cluster_num ) :: distsq
  integer i
  integer it
  integer it_max
  integer j
  integer list(1)
  real, dimension ( dim_num, point_num ) :: point
  integer swap
!
!  Idiot checks.
!
  if ( cluster_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_03 - Fatal error!'
    write ( *, '(a)' ) '  CLUSTER_NUM < 1.'
    stop
  end if

  if ( dim_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_03 - Fatal error!'
    write ( *, '(a)' ) '  DIM_NUM < 1.'
    stop
  end if

  if ( point_num < 1 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_03 - Fatal error!'
    write ( *, '(a)' ) '  POINT_NUM < 1.'
    stop
  end if

  if ( it_max < 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'KMEANS_03 - Fatal error!'
    write ( *, '(a)' ) '  IT_MAX < 0.'
    stop
  end if
!
!  For each observation, calculate the distance from each cluster
!  center, and assign to the nearest.
!
  do i = 1, point_num

    do j = 1, cluster_num
      cluster_energy(j) = sum ( &
        ( point(1:dim_num,i) - cluster_center(1:dim_num,j) )**2 )
    end do

    list = minloc ( cluster_energy(1:cluster_num) )
    cluster(i) = list(1)

  end do
!
!  Determine the cluster populations.
!
  cluster_population(1:cluster_num) = 0
  do i = 1, point_num
    ci = cluster(i)
    cluster_population(ci) = cluster_population(ci) + 1
  end do
!
!  Calculate the mean and sum of squares for each cluster.
!
  cluster_center(1:dim_num,1:cluster_num) = 0.0E+00

  do i = 1, point_num
    j = cluster(i)
    cluster_center(1:dim_num,j) = cluster_center(1:dim_num,j) &
      + point(1:dim_num,i)
  end do

  do i = 1, cluster_num
    cluster_center(1:dim_num,i) = cluster_center(1:dim_num,i) &
      / real ( cluster_population(i) )
  end do
!
!  Carry out the iteration.
!
  it = 0

  do

    if ( it >= it_max ) then
      exit
    end if

    it = it + 1

    swap = 0

    do i = 1, point_num

      ci = cluster(i)

      if ( cluster_population(ci) <= 1 ) then
        cycle
      end if

      do cj = 1, cluster_num

        if ( cj == ci ) then

          distsq(cj) = sum ( &
            ( point(1:dim_num,i) - cluster_center(1:dim_num,cj) )**2 ) &
            * real ( cluster_population(cj) ) &
            / real ( cluster_population(cj) - 1 )

        else if ( cluster_population(cj) == 0 ) then

          cluster_center(1:dim_num,cj) = point(1:dim_num,i)
          distsq(cj) = 0.0E+00

        else

          distsq(cj) = sum ( &
            ( point(1:dim_num,i) - cluster_center(1:dim_num,cj) )**2 ) &
            * real ( cluster_population(cj) ) &
            / real ( cluster_population(cj) + 1 )

        end if

      end do
!
!  Find the index of the minimum value of DISTSQ.
!
      list = minloc ( distsq(1:cluster_num) )
!
!  If that's not the cluster to which point I now belongs, move it there.
!
      if ( list(1) == ci ) then
        cycle
      end if

      cj = list(1)

      cluster_center(1:dim_num,ci) = &
        ( real ( cluster_population(ci) ) * cluster_center(1:dim_num,ci) &
        - point(1:dim_num,i) ) / real ( cluster_population(ci) - 1 )

      cluster_center(1:dim_num,cj) = &
        ( real ( cluster_population(cj) ) * cluster_center(1:dim_num,cj) &
        + point(1:dim_num,i) ) / real ( cluster_population(cj) + 1 )

      cluster_population(ci) = cluster_population(ci) - 1
      cluster_population(cj) = cluster_population(cj) + 1

      cluster(i) = cj

      swap = swap + 1

    end do
!
!  Exit if  no reassignments were made during this iteration.
!
    if ( swap == 0 ) then
      exit
    end if

  end do
!
!  Determine the cluster energies.
!
  call cluster_energy_compute ( dim_num, point_num, cluster_num, point, &
    cluster, cluster_center, cluster_energy )

  return
end
subroutine points_count ( file_in_name, dim_num, point_num )
!
!*******************************************************************************
!
!! POINTS_COUNT counts the valid point coordinates in a file.
!
!
!  Discussion:
!
!    The routine reads every line, and expects to find DIM_NUM
!    real numbers on the line.
!
!    It does not count lines that begin with a comment symbol '#'.
!
!  Modified:
!
!    05 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_IN_NAME, the name of the input file.
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Output, integer POINT_NUM, the number of point coordinate records.
!
  implicit none
!
  integer dim_num
!
  integer bad_num
  integer comment_num
  character ( len = * ) file_in_name
  integer file_in_unit
  integer i
  integer ierror
  integer ios
  character ( len = 100 ) line
  integer point_num
  integer record_num
  real x(dim_num)
!
  call get_unit ( file_in_unit )

  open ( unit = file_in_unit, file = file_in_name, status = 'old', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POINTS_COUNT - Fatal error!'
    write ( *, '(a)' ) '  Could not open the input file: ' // &
      trim ( file_in_name )
    stop
  end if

  comment_num = 0
  point_num = 0
  record_num = 0
  bad_num = 0

  do

    read ( file_in_unit, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      ierror = record_num
      exit
    end if

    record_num = record_num + 1

    if ( line(1:1) == '#' ) then
      comment_num = comment_num + 1
      cycle
    end if

    call s_to_rvec ( line, dim_num, x, ierror )

    if ( ierror == 0 ) then
      point_num = point_num + 1
    else
      bad_num = bad_num + 1
    end if

  end do

  close ( unit = file_in_unit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POINTS_COUNT:'
  write ( *, '(a,i6)' ) '  Number of records:         ', record_num
  write ( *, '(a,i6)' ) '  Number of point records:   ', point_num
  write ( *, '(a,i6)' ) '  Number of comment records: ', comment_num
  write ( *, '(a,i6)' ) '  Number of bad records:     ', bad_num

  return
end
subroutine points_read ( file_in_name, dim_num, point_num, coord )
!
!*******************************************************************************
!
!! POINTS_READ reads point coordinates from a file.
!
!
!  Modified:
!
!    05 March 2002
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) FILE_IN_NAME, the name of the input file.
!
!    Input, integer DIM_NUM, the number of spatial dimensions.
!
!    Input, integer POINT_NUM, the number of points.  The program
!    will stop reading data once POINT_NUM values have been read.
!
!    Output, real COORD(DIM_NUM,POINT_NUM), the point coordinates.
!
  implicit none
!
  integer dim_num
  integer point_num
!
  real coord(dim_num,point_num)
  character ( len = * ) file_in_name
  integer file_in_unit
  integer i
  integer ierror
  integer ios
  character ( len = 100 ) line
  real x(dim_num)
!
  call get_unit ( file_in_unit )

  open ( unit = file_in_unit, file = file_in_name, status = 'old', &
    iostat = ios )

  if ( ios /= 0 ) then
    ierror = 1
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'POINTS_READ - Fatal error!'
    write ( *, '(a)' ) '  Could not open the input file: ' // &
      trim ( file_in_name )
    stop
  end if

  i = 0

  do

    read ( file_in_unit, '(a)', iostat = ios ) line

    if ( ios /= 0 ) then
      ierror = i
      exit
    end if

    if ( line(1:1) == '#' ) then
      cycle
    end if

    call s_to_rvec ( line, dim_num, x, ierror )

    if ( ierror == 0 ) then

      i = i + 1
      if ( i > point_num ) then
        exit
      end if

      coord(1:dim_num,i) = x(1:dim_num)

    end if

  end do

  close ( unit = file_in_unit )

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'POINTS_READ:'
  write ( *, '(a,i6)' ) '  Read coordinate data from file.'

  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:
!
!    19 December 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,
!    and SEED is not changed on output.
!
  implicit none
!
  integer count
  integer count_max
  integer count_rate
  logical, parameter :: debug = .false.
  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

    if ( debug ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANDOM_INITIALIZE'
      write ( *, '(a,i20)' ) '  Initialize RANDOM_NUMBER, user SEED = ', seed
    end if

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    if ( debug ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'RANDOM_INITIALIZE'
      write ( *, '(a,i20)' ) '  Initialize RANDOM_NUMBER, arbitrary SEED = ', &
        seed
    end if

  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 s_to_r ( s, r, ierror, lchar )
!
!*******************************************************************************
!
!! S_TO_R reads a real number from a string.
!
!
!  Discussion:
!
!    This routine will read as many characters as possible until it reaches
!    the end of the string, or encounters a character which cannot be
!    part of the real number.
!
!    Legal input is:
!
!       1 blanks,
!       2 '+' or '-' sign,
!       2.5 spaces
!       3 integer part,
!       4 decimal point,
!       5 fraction part,
!       6 'E' or 'e' or 'D' or 'd', exponent marker,
!       7 exponent sign,
!       8 exponent integer part,
!       9 exponent decimal point,
!      10 exponent fraction part,
!      11 blanks,
!      12 final comma or semicolon.
!
!    with most quantities optional.
!
!  Examples:
!
!    S                 R
!
!    '1'               1.0
!    '     1   '       1.0
!    '1A'              1.0
!    '12,34,56'        12.0
!    '  34 7'          34.0
!    '-1E2ABCD'        -100.0
!    '-1X2ABCD'        -1.0
!    ' 2E-1'           0.2
!    '23.45'           23.45
!    '-4.2E+2'         -420.0
!    '17d2'            1700.0
!    '-14e-2'         -0.14
!    'e2'              100.0
!    '-12.73e-9.23'   -12.73 * 10.0**(-9.23)
!
!  Modified:
!
!    12 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string containing the
!    data to be read.  Reading will begin at position 1 and
!    terminate at the end of the string, or when no more
!    characters can be read to form a legal real.  Blanks,
!    commas, or other nonnumeric data will, in particular,
!    cause the conversion to halt.
!
!    Output, real R, the real value that was read from the string.
!
!    Output, integer IERROR, error flag.
!
!    0, no errors occurred.
!
!    1, 2, 6 or 7, the input number was garbled.  The
!    value of IERROR is the last type of input successfully
!    read.  For instance, 1 means initial blanks, 2 means
!    a plus or minus sign, and so on.
!
!    Output, integer LCHAR, the number of characters read from
!    the string to form the number, including any terminating
!    characters such as a trailing comma or blanks.
!
  implicit none
!
  logical ch_eqi
  character c
  integer ierror
  integer ihave
  integer isgn
  integer iterm
  integer jbot
  integer jsgn
  integer jtop
  integer lchar
  integer nchar
  integer ndig
  real r
  real rbot
  real rexp
  real rtop
  character ( len = * ) s
  character, parameter :: TAB = char ( 9 )
!
  nchar = len_trim ( s )
  ierror = 0
  r = 0.0E+00
  lchar = - 1
  isgn = 1
  rtop = 0.0E+00
  rbot = 1.0E+00
  jsgn = 1
  jtop = 0
  jbot = 1
  ihave = 1
  iterm = 0

  do

    lchar = lchar + 1
    c = s(lchar+1:lchar+1)
!
!  Blank or TAB character.
!
    if ( c == ' ' .or. c == TAB ) then

      if ( ihave == 2 ) then

      else if ( ihave == 6 .or. ihave == 7 ) then
        iterm = 1
      else if ( ihave > 1 ) then
        ihave = 11
      end if
!
!  Comma.
!
    else if ( c == ',' .or. c == ';' ) then

      if ( ihave /= 1 ) then
        iterm = 1
        ihave = 12
        lchar = lchar + 1
      end if
!
!  Minus sign.
!
    else if ( c == '-' ) then

      if ( ihave == 1 ) then
        ihave = 2
        isgn = - 1
      else if ( ihave == 6 ) then
        ihave = 7
        jsgn = - 1
      else
        iterm = 1
      end if
!
!  Plus sign.
!
    else if ( c == '+' ) then

      if ( ihave == 1 ) then
        ihave = 2
      else if ( ihave == 6 ) then
        ihave = 7
      else
        iterm = 1
      end if
!
!  Decimal point.
!
    else if ( c == '.' ) then

      if ( ihave < 4 ) then
        ihave = 4
      else if ( ihave >= 6 .and. ihave <= 8 ) then
        ihave = 9
      else
        iterm = 1
      end if
!
!  Exponent marker.
!
    else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then

      if ( ihave < 6 ) then
        ihave = 6
      else
        iterm = 1
      end if
!
!  Digit.
!
    else if ( ihave < 11 .and. lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      if ( ihave <= 2 ) then
        ihave = 3
      else if ( ihave == 4 ) then
        ihave = 5
      else if ( ihave == 6 .or. ihave == 7 ) then
        ihave = 8
      else if ( ihave == 9 ) then
        ihave = 10
      end if

      call ch_to_digit ( c, ndig )

      if ( ihave == 3 ) then
        rtop = 10.0E+00 * rtop + real ( ndig )
      else if ( ihave == 5 ) then
        rtop = 10.0E+00 * rtop + real ( ndig )
        rbot = 10.0E+00 * rbot
      else if ( ihave == 8 ) then
        jtop = 10 * jtop + ndig
      else if ( ihave == 10 ) then
        jtop = 10 * jtop + ndig
        jbot = 10 * jbot
      end if
!
!  Anything else is regarded as a terminator.
!
    else
      iterm = 1
    end if
!
!  If we haven't seen a terminator, and we haven't examined the
!  entire string, go get the next character.
!
    if ( iterm == 1 .or. lchar+1 >= nchar ) then
      exit
    end if

  end do
!
!  If we haven't seen a terminator, and we have examined the
!  entire string, then we're done, and LCHAR is equal to NCHAR.
!
  if ( iterm /= 1 .and. lchar+1 == nchar ) then
    lchar = nchar
  end if
!
!  Number seems to have terminated.  Have we got a legal number?
!  Not if we terminated in states 1, 2, 6 or 7!
!
  if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then

    ierror = ihave

    return
  end if
!
!  Number seems OK.  Form it.
!
  if ( jtop == 0 ) then
    rexp = 1.0E+00
  else

    if ( jbot == 1 ) then
      rexp = 10.0E+00**( jsgn * jtop )
    else
      rexp = jsgn * jtop
      rexp = rexp / jbot
      rexp = 10.0E+00**rexp
    end if

  end if

  r = isgn * rexp * rtop / rbot

  return
end
subroutine s_to_rvec ( s, n, rvec, ierror )
!
!*******************************************************************************
!
!! S_TO_RVEC reads a real vector from a string.
!
!
!  Modified:
!
!    19 February 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string to be read.
!
!    Input, integer N, the number of values expected.
!
!    Output, real RVEC(N), the values read from the string.
!
!    Output, integer IERROR, error flag.
!    0, no errors occurred.
!    -K, could not read data for entries -K through N.
!
  implicit none
!
  integer n
!
  integer i
  integer ierror
  integer ilo
  integer lchar
  real rvec(n)
  character ( len = * ) s
!
  i = 0

  ilo = 1

  do while ( i < n )

    i = i + 1

    call s_to_r ( s(ilo:), rvec(i), ierror, lchar )

    if ( ierror /= 0 ) then
      ierror = -i
      exit
    end if

    ilo = ilo + lchar

  end do

  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
subroutine word_count ( s, nword )
!
!*******************************************************************************
!
!! WORD_COUNT counts the number of "words" in a string.
!
!
!  Modified:
!
!    14 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, character ( len = * ) S, the string to be examined.
!
!    Output, integer NWORD, the number of "words" in the string.
!    Words are presumed to be separated by one or more blanks.
!
  implicit none
!
  logical blank
  integer i
  integer lens
  integer nword
  character ( len = * ) s
!
  nword = 0
  lens = len ( s )

  if ( lens <= 0 ) then
    return
  end if

  blank = .true.

  do i = 1, lens

    if ( s(i:i) == ' ' ) then
      blank = .true.
    else if ( blank ) then
      nword = nword + 1
      blank = .false.
    end if

  end do

  return
end