if ( 0 < fixed_num ) then    allocate ( fixed(1:dim_num,1:fixed_num) )  else  end if  allocate ( nearest(1:sample_num) )  allocate ( point(1:dim_num,1:n) )  allocate ( point_new(1:dim_num,1:n) )  allocate ( sample(1:dim_num,1:sample_num) )!!  Determine the amount by which the region should be expanded.!  call p00_box ( test, lo, hi )  box_volume = product ( hi(1:3) - lo(1:3) )  h = ( box_volume / real ( n, kind = 8 ) )**( 1.0D+00 / 3.0D+00 )  write ( *, '(a)' ) ' '  write ( *, '(a,i12)' ) '  Number of points requested = ', n  write ( *, '(a,i12)' ) '  Number of CVT iterations = ', iteration_max  write ( *, '(a,i12)' ) '  Number of CVT sample points = ', sample_num  write ( *, '(a,g14.6)' ) '  Using expansion increment H = ', h!!  Get the fixed points.!  if ( 0 < fixed_num ) then    call p00_fixed_points ( test, fixed_num, fixed(1:dim_num,1:fixed_num) )  end if!!  Initialize the sampling points.!  if ( 0 < fixed_num ) then    point(1:dim_num,1:fixed_num) = fixed(1:dim_num,1:fixed_num)  end if  call p00_sample ( test, n-fixed_num, seed, point(1:dim_num,fixed_num+1:n) )  call r8mat_transpose_print_some ( dim_num, n, point, &    1, 1, 2, 10, '  Initial points (first 10 only)' )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  Estimated Voronoi energy:'  write ( *, '(a)' ) ' '!!  Carry out an iteration!  do iteration = 1, iteration_max!!  Get sample points from a "slightly enlarged" region.!    write ( *, * ) 'DEBUG1'    call p00_sample_h1 ( test, sample_num, h, seed, sample )!!  Find the nearest cell generator.!    write ( *, * ) 'DEBUG2'    call find_closest ( dim_num, n, sample_num, point, sample, nearest )!!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).!    write ( *, * ) 'DEBUG3'    point_new(1:dim_num,1:n) = 0.0D+00    count(1:n) = 0    energy = 0.0D+00    do j = 1, sample_num      energy = energy &        + sum ( ( point(1:dim_num,nearest(j)) - sample(1:dim_num,j) )**2 )      point_new(1:dim_num,nearest(j)) = point_new(1:dim_num,nearest(j)) &        + sample(1:dim_num,j)      count(nearest(j)) = count(nearest(j)) + 1    end do!!  Compute the new generators.!  But the fixed points don't move.!    write ( *, * ) 'DEBUG4'    do j = fixed_num+1, n      if ( count(j) /= 0 ) then        point_new(1:dim_num,j) = point_new(1:dim_num,j) &          / real ( count(j), kind = 8 )      end if    end do!!  Project generators back into region.!    write ( *, * ) 'DEBUG5'    call p00_boundary_project ( test, n, point_new )    energy = energy / real ( sample_num, kind = 8 )    write ( *, '(2x,i6,2x,g14.6)' ) iteration, energy!!  Update.!    point(1:dim_num,fixed_num+1:n) = point_new(1:dim_num,fixed_num+1:n)  end do  if ( test < 10 ) then    write ( file_txt_name, '(a,i1,a)' ) 'cvt_p0', test, '_test03.txt'  else    write ( file_txt_name, '(a,i2,a)' ) 'cvt_p', test, '_test03.txt'  end if  call cvt_write ( dim_num, n, seed_start, seed, sample_num, iteration_max, &    energy, point, file_txt_name, comment )  write ( *, '(a)' ) ' '  write ( *, '(a)' ) '  TEST03: wrote data to "' &    // trim ( file_txt_name ) // '".'  deallocate ( count )  if ( 0 < fixed_num ) then    deallocate ( fixed )  end if  deallocate ( nearest )  deallocate ( point )  deallocate ( point_new )   deallocate ( sample )  returnendsubroutine cvt_write ( dim_num, n, seed_start, seed, sample_num, &  iteration_max, energy, point, file_out_name, comment )!*****************************************************************************80!!! CVT_WRITE writes a CVT dataset to a file.!!  Discussion:!!    The initial lines of the file are comments, which begin with a!    '#' character.!!    Thereafter, each line of the file contains the M-dimensional!    components of a CVT generator.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    02 March 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Input, integer DIM_NUM, the spatial dimension.!!    Input, integer N, the number of points.!!    Input, integer SEED_START, the initial random number seed.!!    Input, integer SEED, the current random number seed.!!    Input, integer SAMPLE_NUM, the number of sampling points used on!    each CVT iteration.!!    Input, integer STEP, the number of iterations used in the CVT!    calculation.!!    Input, real ( kind = 8 ) ENERGY, the estimated Voronoi energy.!!    Input, real ( kind = 8 ) CELL_GENERATOR(DIM_NUM,N), the points.!!    Input, character ( len = * ) FILE_OUT_NAME, the name of!    the output file.!!    Input, logical COMMENT, is true if comments may be included.!  implicit none  integer dim_num  integer n  logical comment  real ( kind = 8 ) energy  character ( len = * ) file_out_name  integer file_out_unit  integer ios  integer iteration_max  integer j  real ( kind = 8 ) point(dim_num,n)  integer sample_num  integer seed  integer seed_start  character ( len = 40 ) string  call get_unit ( file_out_unit )  open ( unit = file_out_unit, file = file_out_name, &    status = 'replace', iostat = ios )  if ( ios /= 0 ) then    write ( *, '(a)' ) ' '    write ( *, '(a)' ) 'CVT_WRITE - Fatal error!'    write ( *, '(a)' ) '  Could not open the output file:'    write ( *, '(a)' ) '  "' // trim ( file_out_name ) //'".'    stop  end if  if ( comment ) then    call timestring ( string )    write ( file_out_unit, '(a)'       ) '#  ' // trim ( file_out_name )    write ( file_out_unit, '(a)'       ) &      '#  created by routine CVT_WRITE in CVT_TETRA.F90'    write ( file_out_unit, '(a)'       ) '#  at ' // trim ( string )    write ( file_out_unit, '(a)'       ) '#'    write ( file_out_unit, '(a,i12)'   ) &      '#  Spatial dimension DIM_NUM = ', dim_num    write ( file_out_unit, '(a,i12)'   ) '#  Number of points N =  ', n    write ( file_out_unit, '(a,g14.6)' ) '#  EPSILON (unit roundoff) = ', &      epsilon ( point(1,1) )    write ( file_out_unit, '(a,i12)'   ) '#  Initial SEED = ', seed_start    write ( file_out_unit, '(a,i12)'   ) '#  Current SEED = ', seed    write ( file_out_unit, '(a,i12)'   ) &      '#  Number of sample points =       ', sample_num    write ( file_out_unit, '(a,i12)'   ) &      '#  Number of sampling iterations = ', iteration_max    write ( file_out_unit, '(a,g14.6)' ) &      '#  Estimated Voronoi energy =      ', energy    write ( file_out_unit, '(a)'       )  '#'  end if  write ( string, '(a,i3,a)' ) '(', dim_num, '(f12.6,2x))'  do j = 1, n    write ( file_out_unit, string ) point(1:dim_num,j)  end do  close ( unit = file_out_unit )  returnendsubroutine find_closest ( dim_num, n, sample_num, point, sample, nearest )!*****************************************************************************80!!! FIND_CLOSEST finds the closest point to each sample.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    16 June 2004!!  Author:!!    John Burkardt!!  Parameters:!!    Input, integer DIM_NUM, the spatial dimension.!!    Input, integer N, the number of points.!!    Input, integer SAMPLE_NUM, the number of samples.!!    Input, real ( kind = 8 ) POINT(DIM_NUM,N), the point coordinates.!!    Input, real ( kind = 8 ) SAMPLE(DIM_NUM,SAMPLE_NUM), the sample!    coordinates.!!    Output, integer NEAREST(SAMPLE_NUM), the index of the nearest!    point to each sample.!  implicit none  integer dim_num  integer n  integer sample_num  real ( kind = 8 ) dist  real ( kind = 8 ) dist_min  integer i  integer j  integer nearest(sample_num)  real ( kind = 8 ) point(dim_num,n)  real ( kind = 8 ) sample(dim_num,sample_num)  do i = 1, sample_num    dist_min = huge ( dist_min )    nearest(i) = -1    do j = 1, n      dist = sum ( ( point(1:dim_num,j) - sample(1:dim_num,i) )**2 )      if ( dist < dist_min ) then        dist_min = dist        nearest(i) = j      end if    end do  end do  returnendsubroutine get_unit ( iunit )!*****************************************************************************80!!! 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.!!    If IUNIT = 0, then no free FORTRAN unit could be found, although!    all 99 units were checked (except for units 5, 6 and 9, which!    are commonly reserved for console I/O).!!    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.!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    18 September 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Output, integer IUNIT, the free unit number.!  implicit none  integer i  integer ios  integer iunit  logical lopen  iunit = 0  do i = 1, 99    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) 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  returnendsubroutine timestring ( string )!*****************************************************************************80!!! TIMESTRING writes the current YMDHMS date into a string.!!  Example:!!    STRING = '31 May 2001   9:45:54.872 AM'!!  Licensing:!!    This code is distributed under the GNU LGPL license. !!  Modified:!!    06 August 2005!!  Author:!!    John Burkardt!!  Parameters:!!    Output, character ( len = * ) STRING, contains the date information.!    A character length of 40 should always be sufficient.!  implicit none  character ( len = 8 ) ampm  integer d  integer h  integer m  integer mm  character ( len = 9 ), parameter, dimension(12) :: month = (/ &    'January  ', 'February ', 'March    ', 'April    ', &    'May      ', 'June     ', 'July     ', 'August   ', &    'September', 'October  ', 'November ', 'December ' /)  integer n  integer s  character ( len = * ) string  integer values(8)  integer y  call date_and_time ( values = values )  y = values(1)  m = values(2)  d = values(3)  h = values(5)  n = values(6)  s = values(7)  mm = values(8)  if ( h < 12 ) then    ampm = 'AM'  else if ( h == 12 ) then    if ( n == 0 .and. s == 0 ) then      ampm = 'Noon'    else      ampm = 'PM'    end if  else    h = h - 12    if ( h < 12 ) then      ampm = 'PM'    else if ( h == 12 ) then      if ( n == 0 .and. s == 0 ) then        ampm = 'Midnight'      else        ampm = 'AM'      end if    end if  end if  write ( string, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )  returnend