subroutine cvt_size_iteration ( dim_num, box_min, box_max, cell_num, &
  cell_generator, cell_weight, sample_num, cell_volume, region_volume )
!
!******************************************************************************
!
!! CVT_SIZE_ITERATION takes one step of the CVT size iteration.
!
!
!  Discussion:
!
!    The routine is given a set of points, called "generators", which
!    define a tessellation of the region into Voronoi cells.  Each point
!    defines a cell.  Each cell, in turn, has a centroid, but it is
!    unlikely that the centroid and the generator coincide.
!
!    Each time this CVT iteration is carried out, an attempt is made
!    to modify the generators in such a way that they are closer and
!    closer to being the centroids of the Voronoi cells they generate.
!
!    A large number of sample points are generated, and the nearest generator
!    is determined.  A count is kept of how many points were nearest to each
!    generator.  Once the sampling is completed, the location of all the
!    generators is adjusted.  This step should decrease the discrepancy
!    between the generators and the centroids.
!
!  Modified:
!
!    16 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, real CELL_WEIGHT(CELL_NUM), the cell weights.
!
!    Input, integer SAMPLE_NUM, the number of sample points.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(1:dim_num)
  real box_min(1:dim_num)
  real c1
  real c2
  integer cell_count(cell_num)
  real cell_generator(dim_num,cell_num)
  real cell_generator2(dim_num,cell_num)
  real cell_volume(cell_num)
  real cell_weight(cell_num)
  logical, parameter :: debug = .false.
  integer i
  integer j
  integer k
  integer nearest
  integer ngen
  real region_volume
  integer sample_num
  real x(dim_num)
!
  cell_generator2(1:dim_num,1:cell_num) = 0.0E+00
  cell_count(1:cell_num) = 0

  do j = 1, sample_num
!
!  Generate a sampling point X.
!
    call region_sampler ( dim_num, box_min, box_max, x, ngen )
!
!  Find the nearest cell generator.
!
    call find_closest_size ( dim_num, x, cell_num, cell_generator, &
      cell_weight, nearest )
!
!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).
!
    if ( nearest < 1 .or. nearest > cell_num ) then
      write ( *, * ) ' '
      write ( *, * ) 'NEAREST = ',nearest
      stop
    end if

    cell_generator2(1:dim_num,nearest) = &
      cell_generator2(1:dim_num,nearest) + x(1:dim_num)

    cell_count(nearest) = cell_count(nearest) + 1

  end do
!
!  Compute the new generators.
!
  do j = 1, cell_num

    if ( cell_count(j) /= 0 ) then

      cell_generator(1:dim_num,j) = cell_generator2(1:dim_num,j) &
        / real ( cell_count(j) )

    end if

  end do

  cell_volume(1:cell_num) = cell_count(1:cell_num) * region_volume &
    / real ( sample_num )

  return
end
subroutine find_closest_size ( dim_num, x, cell_num, cell_generator, &
  cell_weight, nearest )
!
!******************************************************************************
!
!! FIND_CLOSEST_SIZE 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:
!
!    03 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, real X(DIM_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.
!
!    Input, real CELL_WEIGHT(CELL_NUM), the cell weights.
!
!    Output, integer NEAREST, the index of the nearest cell generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real cell_generator(dim_num,cell_num)
  real cell_weight(cell_num)
  real distance
  real dist_sq
  integer i
  integer nearest
  real x(dim_num)
!
  nearest = 0
  distance = huge ( distance )

  do i = 1, cell_num

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

    dist_sq = sqrt ( dist_sq ) / cell_weight(i)

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

  end do

  return
end
subroutine generator_init ( dim_num, box_min, box_max, cell_num, cell_generator )
!
!******************************************************************************
!
!! GENERATOR_INIT initializes the Voronoi cell generators.
!
!
!  Discussion:
!
!    The points initialized here will be used to generate a tessellation
!    of the region into Voronoi cells.  Each generator point defines a
!    cell.  The CVT algorithm will try to modify these initial generators
!    in such a way that they are also the centroids of the cells they generate.
!
!    It is probably better to use Halton points for the centroids than
!    uniform random values, in the sense that the algorithm will probably
!    converge more quickly.
!
!  Modified:
!
!    03 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.
!
!    Output, real CELL_GENERATOR(DIM_NUM,CELL_NUM), the Voronoi cell
!    generators.
!
  implicit none
!
  integer cell_num
  integer dim_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  real cell_generator(dim_num,cell_num)
  integer i
  integer ngen
!
  do i = 1, cell_num

    call region_sampler ( dim_num, box_min, box_max, cell_generator(1,i), ngen )

  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, box_min, box_max, x, 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, either by a uniform random
!    number generator, or from a vector Halton sequence.
!
!    If a user-supplied routine determines that this point is
!    within the physical region, this routine returns.  Otherwise,
!    a new random point is chosen.
!
!    The entries of the local vector HALTON_BASE should be distinct primes.
!    Right now, we're assuming DIM_NUM is no greater than 3.
!
!  Modified:
!
!    03 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.
!
!    Output, real X(DIM_NUM), the random point.
!
!    Output, integer NGEN, the number of points that were generated.
!    This is at least 1, but may be larger if some points were rejected.
!
  integer dim_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  integer i
  integer ival
  integer ngen
  real r(dim_num)
  real x(dim_num)
!
  ngen = 0

  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)' ) x(1:dim_num)
      stop
    end if

    call random_number ( r(1:dim_num) )
!
!  Determine a point in the bounding box.
!
    x(1:dim_num) = ( ( 1.0E+00 - r(1:dim_num) ) * box_min(1:dim_num) &
                            + r(1:dim_num)   * box_max(1:dim_num) )

    call test_region ( x, dim_num, ival )

    if ( ival == 1 ) then
      exit
    end if

  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