subroutine cell_volume_computation ( dim_num, box_min, box_max, cell_num, &  
  cell_generator, sample_num, cell_volume, region_volume )
!
!******************************************************************************
!
!! CELL_VOLUME_COMPUTATION estimates the cell volumes by sampling.
!
!
!  Modified:
!
!    25 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer CELL_NUM, the number of Voronoi cells.
!
!    Input/output, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the Voronoi
!    cell generators.  On output, these have been modified.
!
  implicit none
!
  integer cell_num
  integer dim_num
  integer sample_num
!
  real box_max(1:dim_num)
  real box_min(1:dim_num)
  integer cell_count(cell_num)
  real cell_generator(dim_num,cell_num)
  real cell_volume(cell_num)
  integer i
  integer j
  integer nearest(sample_num)
  integer ngen
  real region_volume
  real sample(dim_num,sample_num)
!
!  Generate the sampling points.
!
  call region_sampler ( dim_num, sample_num, box_min, box_max, sample, ngen )
!
!  Find closest generators.
!
  call find_closest ( dim_num, sample_num, sample, cell_num, &
    cell_generator, nearest )

  cell_count(1:cell_num) = 0
  do i = 1, sample_num
    j = nearest(i)
    if ( j < 1 .or. j > cell_num ) then
      write ( *, * ) 'i = ', i
      write ( *, * ) 'j = ', j
      write ( *, * ) 'Cell_num = ', cell_num
      stop
    end if
    cell_count(j) = cell_count(j) + 1
  end do

  do j = 1, cell_num
    cell_volume(j) = region_volume * real ( cell_count(j) ) &
      / real ( sample_num )
  end do

  return
end
subroutine cvt_max_iteration ( dim_num, box_min, box_max, cell_num, &  
  cell_generator, cell_count_desired, sample_num, cell_volume, region_volume )
!
!******************************************************************************
!
!! CVT_MAX_ITERATION takes one step of the CVT Max iteration.
!
!
!  Discussion:
!
!    On input, each generator has been assigned a desired count
!    or capacity.  The sum of these capacitities determines the
!    number of sample points used.  These sample points are then
!    assigned to the generators according to the following rules:
!
!    * Each generator should get exactly its capacity of points.
!
!    * the points which are closest to a generator G are assigned to
!      that generator, in order of closeness, up to the capacity of
!      the generator.  These are the initial assignments.
!
!    * Once initial assignments have been made, there will generally
!      be a number of unassigned sample points remaining.  The assignment
!      process is repeated, but now using only the generators that
!      still have open slots.  Each unassigned sample point is assigned 
!      to its nearest generator, in order of closeness, up to the
!      capacity of that generator.
!
!    * This process is repeated until all points have been assigned.
!
!    Once the sampling and assignment is completed, the location of all 
!    the generators is adjusted.  This step should decrease the discrepancy
!    between the generators and the centroids.
!
!  Modified:
!
!    23 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer CELL_NUM, the number of Voronoi cells.
!
!    Input/output, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the Voronoi
!    cell generators.  On output, these have been modified.
!
!    Input, integer CELL_SIZE(CELL_NUM), the maximum number of points
!    that can be accepted by each cell generator.  The total of these
!    entries is the number of sample points that will be generated.
!
  implicit none
!
  integer cell_num
  integer dim_num
  integer sample_num
!
  real box_max(1:dim_num)
  real box_min(1:dim_num)
  integer candidate_indx(sample_num)
  integer candidate_num
  integer cell_count(cell_num)
  integer cell_count_desired(cell_num)
  real cell_generator(dim_num,cell_num)
  real cell_volume(cell_num)
  logical, parameter :: debug = .true.
  real distance(sample_num)
  integer i
  integer indx(sample_num)
  integer it
  integer, parameter :: it_max = 3
  integer jj
  integer j
  integer k
  integer nearest(sample_num)
  integer ngen
  integer opening_num
  real region_volume
  real sample(dim_num,sample_num)
  integer sample_cell(sample_num)
!
  cell_count(1:cell_num) = 0
  sample_cell(1:sample_num) = 0
  nearest(1:sample_num) = 0
  distance(1:sample_num) = huge ( 1.0E+00 )
!
!  Generate the sampling points.
!
  call region_sampler ( dim_num, sample_num, box_min, box_max, sample, ngen )
!
!  Then, for each generator, take the candidates in order of
!  closeness, up to the capacity of the generator.
!
  do it = 1, it_max

    if ( debug ) then
      write ( *, * ) ' '
      write ( *, * ) '  MAX-Iteration number ', it
      write ( *, * ) ' '
    end if

    do i = 1, sample_num
      if ( sample_cell(i) == 0 ) then
        nearest(i) = 0
      end if
    end do
!
!  For each unassigned sample point, compute the nearest available generator.
!
    call find_closest_available ( dim_num, sample_num, sample, cell_num, &
      cell_generator, cell_count, cell_count_desired, nearest, distance )

    if ( debug ) then
    write ( *, * ) ' '
    write ( *, * ) '  I, NEAREST, SAMPLE_CELL, DISTANCE'
    write ( *, * ) ' '
    do i = 1, sample_num
      write ( *, '(i5,i5,i5,g14.6)' ) &
        i, nearest(i), sample_cell(i), distance(i)
    end do
    end if
!
!  Now, for each unfilled generator, consider the unclaimed points 
!  that view that generator as their nearest neighbor.  Sort them
!  by distance, and assign as many as possible to the generator.
!  
    do j = 1, cell_num

      opening_num = cell_count_desired(j) - cell_count(j)

      if ( opening_num > 0 ) then

        if ( debug ) then
        write ( *, * ) ' '
        write ( *, * ) 'Generator = ', j
        end if

        candidate_num = 0

        do i = 1, sample_num

          if ( sample_cell(i) == 0 .and. nearest(i) == j ) then
            candidate_num = candidate_num + 1
            candidate_indx(candidate_num) = i
          end if
          
        end do
!
!  Now sort the index entries by distance.
!
        call rvec_sort_heap_mask_a ( sample_num, distance, candidate_num, &
          candidate_indx, indx )

        do k = 1, min ( candidate_num, opening_num )
          sample_cell(candidate_indx(indx(k))) = j
          cell_count(j) = cell_count(j) + 1
        end do
!
!  Let's see what we have.
!
        if ( debug ) then
        write ( *, * ) ' '
        write ( *, * ) '  I, SAMPLE_CELL(I), DISTANCE(I)'
        write ( *, * ) ' '

        do i = 1, sample_num
          write ( *, '(2i6,g14.6)' ) i, sample_cell(i), distance(i)
        end do

        write ( *, * ) ' '
        write ( *, * ) '  J, CELL_COUNT(J), CELL_COUNT_DESIRED(J)'
        write ( *, * ) ' '
        do jj = 1, cell_num
          write ( *, '(3i6)' ) jj, cell_count(jj), cell_count_desired(jj)
        end do
        end if

      end if

    end do

  end do
!
!  Now compute the new generators.
!
  cell_generator(1:dim_num,1:cell_num) = 0.0

  do i = 1, sample_num
    j = sample_cell(i)
    cell_generator(1:dim_num,j) = cell_generator(1:dim_num,j) + &
      sample(1:dim_num,i)
  end do

  do j = 1, cell_num
    cell_generator(1:dim_num,j) = cell_generator(1:dim_num,j) / real ( cell_count(j) )
  end do

  do j = 1, cell_num
    cell_volume(j) = region_volume * real ( cell_count(j) ) / real ( sample_num )
  end do

  return
end
subroutine find_closest ( dim_num, sample_num, sample, cell_num, cell_generator, &
  nearest )
!
!******************************************************************************
!
!! FIND_CLOSEST finds the Voronoi cell generator closest to a point X.
!
!
!  Discussion:
!
!    This routine finds the closest Voronoi cell generator by checking every
!    one.  For problems with many cells, this process can take the bulk
!    of the CPU time.  Other approaches, which group the cell generators into
!    bins, can run faster by a large factor.
!
!  Modified:
!
!    18 January 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real SAMPLE(DIM_NUM,SAMPLE_NUM), the point to be checked.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Output, integer NEAREST(SAMPLE_NUM), the index of the nearest 
!    cell generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
  integer sample_num
!
  real cell_generator(dim_num,cell_num)
  real distance
  real dist_sq
  integer i
  integer j
  integer nearest(sample_num)
  real sample(dim_num,sample_num)
!
  nearest(1:sample_num) = 0

  do i = 1, sample_num

    distance = huge ( distance )

    do j = 1, cell_num

      dist_sq = sum ( ( cell_generator(1:dim_num,j) - sample(1:dim_num,i) )**2 )

      if ( dist_sq < distance ) then
        distance = dist_sq
        nearest(i) = j
      end if

    end do

    distance = sqrt ( distance )

  end do

  return
end
subroutine find_closest_available ( dim_num, sample_num, sample, cell_num, &
  cell_generator, cell_count, cell_count_desired, nearest, distance )

!
!******************************************************************************
!
!! FIND_CLOSEST_AVAILABLE matches points to nearest available generators.
!
!
!  Discussion:
!
!    Do we pass in the whole array of sample points and generators
!    every time?  If so, how do we mark the NEAREST results as
!    tentative?  Because somebody has to take all the tentative 
!    results, sort them by distance, and take no more than the limit.
!
!  Modified:
!
!    22 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer SAMPLE_NUM, the number of sample points.
!
!    Input, real SAMPLE(DIM_NUM,SAMPLE_NUM), the sample points.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the cell generators.
!
!    Input/output, integer CELL_COUNT(CELL_NUM), the number of points
!    associated with each cell generator.  On return, the entry
!    for the generator associated with the input point X will have
!    been incremented.
!
!    Input, integer CELL_COUNT_DESIRED(CELL_NUM), the maximum number 
!    of points that can be accepted by each cell generator.
!
!    Output, integer NEAREST(SAMPLE_NUM), the index of the nearest 
!    "available" cell generator.  
!
!    Output, real DISTANCE(SAMPLE_NUM), the distance between each sample
!    point and its nearest available cell generator.
!
  implicit none
!
  integer cell_num
  integer dim_num
  integer sample_num
!
  integer cell_count(cell_num)
  integer cell_count_desired(cell_num)
  real cell_generator(dim_num,cell_num)
  real distance(sample_num)
  real dist
  integer i
  integer j
  integer nearest(sample_num)
  real sample(dim_num,sample_num)
!
  do i = 1, sample_num

    if ( nearest(i) == 0 ) then

      distance(i) = huge ( distance(i) )

      do j = 1, cell_num

        if ( cell_count(j) < cell_count_desired(j) ) then

          dist = sqrt ( &
            sum ( ( cell_generator(1:dim_num,j) - sample(1:dim_num,i) )**2 ) )

          if ( dist < distance(i) ) then 
            distance(i) = dist
            nearest(i) = j
          end if

        end if

      end do

    end if

  end do

  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.
!
  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 ( *, * ) ' '
    write ( *, * ) 'RANDOM_INITIALIZE'
    write ( *, * ) '  Initialize RANDOM_NUMBER with user SEED = ', seed

  else

    call system_clock ( count, count_rate, count_max )

    seed = count

    write ( *, * ) ' '
    write ( *, * ) 'RANDOM_INITIALIZE'
    write ( *, * ) '  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 region_sampler ( dim_num, sample_num, box_min, box_max, sample, ngen )
!
!******************************************************************************
!
!! REGION_SAMPLER returns a sample point in the physical region.
!
!
!  Discussion:
!
!    The calculations are done in DIM_NUM dimensional space.
!
!    The physical region is enclosed in a bounding box.
!
!    A point is chosen in the bounding box by a uniform random
!    number generator.
!
!    If a user-supplied routine determines that this point is
!    within the physical region, this routine returns.  Otherwise,
!    a new random point is chosen.
!
!  Modified:
!
!    22 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer SAMPLE_NUM, the number of points to be generated.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Output, real SAMPLE(DIM_NUM,SAMPLE_NUM), the sample points.
!
!    Output, integer NGEN, the number of points that were generated.
!    This is at least SAMPLE_NUM, but may be larger if some points 
!    were rejected.
!
  integer dim_num
  integer sample_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  integer ival
  integer j
  integer ngen
  real r(dim_num)
  real sample(dim_num,sample_num)
!
  ngen = 0

  do j = 1, sample_num

    do

      ngen = ngen + 1

      if ( ngen > 10000 ) then
        write ( *, * ) ' '
        write ( *, * ) 'REGION_SAMPLER - Fatal error!'
        write ( *, * ) '  Generated ', ngen, ' rejected points in a row.'
        write ( *, * ) '  There may be a problem with the geometry definition.'
        write ( *, * ) ' '
        write ( *, * ) '  Current random value is:'
        write ( *, '(3g14.6)' ) r(1:dim_num)
        write ( *, * ) ' '
        write ( *, * ) '  Current sample point is:'
        write ( *, '(3g14.6)' ) sample(1:dim_num,j)
        stop
      end if
!
!  Generate a point at random.
!
      call random_number ( r(1:dim_num) )
!
!  Determine a point in the bounding box.
!
      sample(1:dim_num,j) = ( ( 1.0E+00 - r(1:dim_num) ) * box_min(1:dim_num) &
                                     + r(1:dim_num)   * box_max(1:dim_num) )

      call test_region ( sample(1:dim_num,j), dim_num, ival )

      if ( ival == 1 ) then
        exit
      end if

    end do

  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
!
  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