program cvt_main
!
!******************************************************************************
!
!! CVT_MAIN runs a test of CVT, the centroidal Voronoi tesselation library.
!
!
!  Variables:
!
!
!    Geometry Variables:
!    ------------------
!
!    Input, integer NDIM, the spatial dimension, which should be 2 or 3.
!
!    double precision BOX_MIN(NDIM), BOX_MAX(NDIM), the minimum and maximum 
!    coordinates of the bounding box.  
!
!
!    CVT Algorithm Variables:
!    -----------------------
!
!    integer N, the number of Voronoi cells to generate.
!    A typical value is 256.
!
!    integer MAXIT, the maximum number of correction iterations used in the
!    Voronoi calculation.  A typical value is 10.
!
!    integer NS_CVT, the average number of sampling points tested per 
!    Voronoi cell, on each step of the correction iteration.  A typical 
!    value is 5000.
!
!    Input, integer RANDOM_GENERATOR, specifies how the Voronoi cell 
!    generators are to be initialized.
!    0, use the F90 RANDOM_NUMBER routine;
!    1, use the Halton sequence. (PREFERRED)
!
!    double precision CELL_GENERATOR(NDIM,N), the Voronoi cell generators 
!    of the Voronoi tesselation, as approximated by the CVT algorithm.  This
!    is the output quantity of most interest.
!
!
!    Moment Calculation Variables:
!    ----------------------------
!
!    integer NS_MOM, the average number of sampling points tested per 
!    Voronoi cell, for the moment calculation.  A typical value is 5000.
!
!    logical REGION_VOLUME_GIVEN, 
!    0, the area or volume is already available in REGION.
!    nonzero, the area or volume needs to be calculated.
!
!    double precision REGION_VOLUME.
!    If REGION_VOLUME_GIVEN, then REGION_VOLUME must be input by the user.
!    Otherwise, REGION_VOLUME is approximated computationally.
!
!    double precision CELL_VOLUME(N), the volume of the Voronoi cells.
!
!    double precision CELL_CENTROID(NDIM,N), the centroids of the Voronoi cells.
!
!    double precision CELL_MOMENT(NDIM,NDIM,N), the second moments of the
!    Voronoi cells.
!
!
!  Miscellaneous Variables:
!  -----------------------
!
!    integer SEED, determines how to initialize the RANDOM_NUMBER routine.
!    If SEED is zero on input, then RANDOM_INITIALIZE will make up a seed
!    from the current real time clock reading.
!    If SEED is nonzero, then a reproducible sequence of random numbers
!    defined by SEED will be chosen.
!
  integer, parameter :: n = 16
  integer, parameter :: ndim = 2
!
  double precision, parameter, dimension ( ndim ) :: box_max = &
    (/ 100.0D+00, 100.0D+00 /)
! double precision, parameter, dimension ( ndim ) :: box_max = &
!   (/ 100.0D+00, 100.0D+00, 20.0D+00 /)
  double precision, parameter, dimension ( ndim ) :: box_min = &
    (/ 0.0D+00, 0.0D+00 /)
! double precision, parameter, dimension ( ndim ) :: box_min = &
!   (/ 0.0D+00, 0.0D+00, 0.0D+00 /)
  double precision cell_centroid(ndim,n)
  double precision cell_generator(ndim,n)
  double precision cell_moment(ndim,ndim,n)
  double precision cell_volume(n)
  integer clock0
  integer clock1
  integer clock_max
  integer clock_rate
  character ( len = 8 ) date
  real etime0
  real etime1
  integer i
  integer ios
  integer it
  integer, parameter :: maxit = 10
  integer, parameter :: ns_cvt = 1000
  integer, parameter :: ns_mom = 1000
  integer, parameter :: random_generator = 1
  double precision, parameter :: region_volume = 1700.0D+00
! double precision, parameter :: region_volume = 20.0D+00 * 1700.0D+00
  logical, parameter :: region_volume_given = .true.
  integer, save :: seed = 0
  real tarray(2)
  character ( len = 10 ) time
  real time0
  real time1
  integer, dimension ( n ) :: updates
!
!  Print introduction and options.
!
  call date_and_time ( date, time )

  write ( *, * ) ' '
  write ( *, * ) 'CVT_MAIN'
  write ( *, * ) '  A sample problem for the probabilistic'
  write ( *, * ) '  Centroidal Voronoi Tesselation algorithm.'
  write ( *, * ) ' '
  write ( *, * ) '  Today''s date: ', date
  write ( *, * ) '  Today''s time: ', time
  write ( *, * ) ' '
  write ( *, * ) '  Given a region in 2D or 3D, the problem is to determine'
  write ( *, * ) '  GENERATORS, a set of points which define a division'
  write ( *, * ) '  of the region into Voronoid cells, which are also CENTROIDS'
  write ( *, * ) '  of the Voronoi cells.'

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) 'Geometry parameters:'
  write ( *, * ) '-------------------'
  write ( *, * ) ' '
  write ( *, * ) '  The spatial dimension is NDIM = ', ndim
  write ( *, * ) ' '
  write ( *, * ) '  The minimum corner of the bounding box is:'
  write ( *, * ) box_min(1:ndim)
  write ( *, * ) '  The maximum corner of the bounding box is:'
  write ( *, * ) box_max(1:ndim)

  write ( *, * ) ' '
  write ( *, * ) 'CVT Algorithm parameters:'
  write ( *, * ) '-------------------------'
  write ( *, * ) ' '
  write ( *, * ) '  The number of Voronoi cells to generate: ', n
  write ( *, * ) '  Number of iterations to determine CVT: ', maxit
  write ( *, * ) '  Number of sampling points per Voronoi cell: ', ns_cvt

  if ( random_generator == 0 ) then
    write ( *, * ) '  Voronoi cell generators are initialized by RANDOM_NUMBER.'
  else if ( random_generator == 1 ) then
    write ( *, * ) '  Voronoi cell generators are initialized by Halton.'
  end if

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) 'Moment parameters:'
  write ( *, * ) '------------------'
  write ( *, * ) ' '
  write ( *, * ) '  Number of sampling points per Voronoi cell: ', ns_mom
  if ( region_volume_given ) then
    write ( *, * ) ' '
    write ( *, * ) '  The volume of the region is given.'
    write ( *, * ) '  It is specified as REGION_VOLUME = ', region_volume
  else
    write ( *, * ) ' '
    write ( *, * ) '  The volume of the region is not given.'
    write ( *, * ) '  It will be estimated.'
  end if
!
!  Initialize the random number generator.
!
!  If SEED is zero on input, then the routine will make up a seed.
!  If SEED is nonzero, then a reproducible sequence of random numbers
!  defined by SEED will be chosen.
!
  call random_initialize ( seed )
!
!  Begin timing.
!
!  The F90 CPU_TIME routine, at least on the DEC_ALPHA, can "wrap around"
!  after a relatively short time, giving negative timings.  So we also
!  call ETIME, which should be measuring the same thing, for checking.
!
!  The F90 SYSTEM_CLOCK routine should measure real execution time.
!
  call system_clock ( clock0, clock_rate, clock_max )
  call cpu_time ( time0 )
  call etime ( tarray )
  etime0 = tarray(1) + tarray(2)
!
!  Initialize the Voronoi cell generators.
!
  call generator_init ( ndim, box_min, box_max, n, cell_generator, &
    random_generator )
!
!  Carry out the CVT iteration, which drives the Voronoi generators 
!  and Voronoi centroids closer and closer.
!
  updates(1:n) = 1

  do it = 1, maxit

    call cvt_iteration ( ndim, box_min, box_max, n, cell_generator, ns_cvt, &
      updates )

  end do
!
!  Calculate cell moments.
!
  call vcm ( ndim, box_min, box_max, n, cell_generator, ns_mom, &
    region_volume_given, region_volume, cell_volume, cell_centroid, cell_moment )
!
!  End timing the "interesting" part.
!
  call etime ( tarray )
  etime1 = tarray(1) + tarray(2)
  call cpu_time ( time1 )
  call system_clock ( clock1, clock_rate, clock_max )

  write ( *, * ) ' '
  write ( *, * ) 'Elapsed CPU time, CPU_TIME: ', time1-time0, ' seconds.'
  write ( *, * ) 'Elapsed CPU time, ETIME:    ', etime1-etime0, ' seconds.'
  write ( *, * ) 'Elapsed time, SYSTEM_CLOCK: ', real ( clock1 - clock0 ) &
    / real ( clock_rate ), ' seconds.'
!
!  Compute quality checks.
!
  call quality ( ndim, n, cell_moment, cell_volume, region_volume )
!
!  Write generators and moments to files.
!
  open ( unit = 1, file = 'cvt_generators.dat', status = 'replace', &
    iostat = ios )
  do i = 1, n
    write ( 1, '(3g15.8)' ) cell_generator(1:ndim,i)
  end do
  close ( unit = 1 )

  open ( unit = 1, file = 'cvt_volume.dat', status = 'replace', iostat = ios )
  write ( 1, '(g15.8)' ) cell_volume(1:n)
  close ( unit = 1 )

  open ( unit = 1, file = 'cvt_centroid.dat', status = 'replace', iostat = ios )
  do i = 1, n
    write ( 1,'(3g15.8)' ) cell_centroid(1:ndim,i)
  end do
  close ( unit = 1 )

  open ( unit = 1, file = 'cvt_moment.dat', status = 'replace', iostat = ios )
  do i = 1, n
    write ( 1, '(3g15.8)' ) cell_moment(1:ndim,1:ndim,i)
  end do
  close ( unit = 1 )

  write ( *, * ) ' '
  write ( *, * ) 'CVT_MAIN'
  write ( *, * ) '  Normal end of execution.'

  stop
end
subroutine cvt_iteration ( ndim, box_min, box_max, n, cell_generator, ns_cvt, &
  updates )
!
!******************************************************************************
!
!! CVT_ITERATION takes one step of the CVT iteration.
!
!
!  Discussion:
!
!    The routine is given an set of points, called "generators", which
!    define a tesselation 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 April 2001
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision BOX_MIN(NDIM), BOX_MAX(NDIM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer N, the number of Voronoi cells.
!
!    Input/output, double precision CELL_GENERATOR(NDIM,N), the Voronoi
!    cell generators.  On output, these have been modified
!
!    Input, integer NS_CVT, the number of sample points per generator.
!
!    Input/output, integer UPDATES(N), counts the number of times a cell
!    center has been updated.  Before the first call, all the entries of
!    UPDATES should be set to 1.  After each iteration, UPDATES will be
!    incremented by 1 for each cell generator that was updated.  Normally,
!    all of them will be so updated.
!
  integer n
  integer ndim
!
  double precision box_max(1:ndim)
  double precision box_min(1:ndim)
  double precision c1
  double precision c2
  double precision cell_generator(ndim,n)
  double precision cell_generator2(ndim,n)
  integer count(n)
  integer i
  integer j
  integer k
  integer nearest
  integer ngen
  integer ns_cvt
  integer, parameter :: random_sampler = 0
  logical reset
  integer updates(n)
  double precision x(ndim)
!
  cell_generator2(1:ndim,1:n) = 0.0D+00
  count(1:n) = 0
  reset = .false.

  do j = 1, ns_cvt * n
!
!  Generate a sampling point X.
!
    call region_sampler ( ndim, box_min, box_max, x, &
     random_sampler, reset, ngen )
!
!  Find the nearest cell generator.
!
    call find_closest ( ndim, x, n, cell_generator, nearest )
!
!  Add X to the averaging data for CELL_GENERATOR(*,NEAREST).
!
    cell_generator2(1:ndim,nearest) = &
      cell_generator2(1:ndim,nearest) + x(1:ndim)

    count(nearest) = count(nearest) + 1

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

    if ( count(j) /= 0 ) then

      cell_generator(1:ndim,j) = cell_generator2(1:ndim,j) / dble ( count(j) )
      updates(j) = updates(j) + 1

    end if

  end do

  return
end
function dmat_det_2d ( a )
!
!*******************************************************************************
!
!! DMAT_DET_2D computes the determinant of a 2 by 2 matrix.
!
!
!  Formula:
!
!    The determinant of a 2 by 2 matrix is
!
!      a11 * a22 - a12 * a21.
!
!  Modified:
!
!    01 March 1999
!
!  Author:
!
!    John Burkardt,
!    Department of Mathematics,
!    Iowa State University.
!
!  Parameters:
!
!    Input, double precision A(2,2), the matrix whose determinant is desired.
!
!    Output, double precision RMAT_DET_2D, the determinant of the matrix.
!
  double precision a(2,2)
  double precision dmat_det_2d
!
  dmat_det_2d = a(1,1) * a(2,2) - a(1,2) * a(2,1)

  return
end
function dmat_det_3d ( a )
!
!*******************************************************************************
!
!! DMAT_DET_3D computes the determinant of a 3 by 3 matrix.
!
!
!  Formula:
!
!    The determinant of a 3 by 3 matrix is
!
!        a11 * a22 * a33 - a11 * a23 * a32
!      + a12 * a23 * a31 - a12 * a21 * a33
!      + a13 * a21 * a32 - a13 * a22 * a31
!
!  Modified:
!
!    01 March 1999
!
!  Author:
!
!    John Burkardt,
!    Department of Mathematics,
!    Iowa State University.
!
!  Parameters:
!
!    Input, double precision A(3,3), the matrix whose determinant is desired.
!
!    Output, double precision RMAT_DET_3D, the determinant of the matrix.
!
  double precision a(3,3)
  double precision dmat_det_3d
!
  dmat_det_3d = &
         a(1,1) * ( a(2,2) * a(3,3) - a(2,3) * a(3,2) ) &
       + a(1,2) * ( a(2,3) * a(3,1) - a(2,1) * a(3,3) ) &
       + a(1,3) * ( a(2,1) * a(3,2) - a(2,2) * a(3,1) )

  return
end
subroutine find_closest ( ndim, x, n, 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
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision X(NDIM), the point to be checked.
!
!    Input, integer N, the number of cell generatorrs.
!
!    Input, double precision CELL_GENERATOR(NDIM,N), the cell generators.
!
!    Output, integer NEAREST, the index of the nearest cell generators.
!
  integer n
  integer ndim
!
  double precision cell_generator(ndim,n)
  double precision distance
  double precision dist_sq
  integer i
  integer nearest
  double precision x(ndim)
!
  nearest = 0
  distance = huge ( distance )

  do i = 1, n

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

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

  end do

  distance = sqrt ( distance )

  return
end
subroutine generator_init ( ndim, box_min, box_max, n, cell_generator, &
  random_generator )
!
!******************************************************************************
!
!! GENERATOR_INIT initializes the Voronoi cell generators.
!
!
!  Discussion:
!
!    The points initialized here will be used to generate a tesselation
!    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:
!
!    02 April 2001
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision BOX_MIN(NDIM), BOX_MAX(NDIM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer N, the number of Voronoi cells.
!
!    Output, double precision CELL_GENERATOR(NDIM,N), the Voronoi cell
!    generators.
!
!    Input, logical USE_DIATOM, is TRUE if DIATOM is to be called
!    to determine whether a point lies in the physical region; if it is
!    FALSE than a much simplified routine is used.
!
!    Input, double precision DR, a tolerance used by DIATOM when testing
!    whether a point is within, outside of, or on the boundary of the
!    physical region.
!
!    Input, integer RANDOM_GENERATOR, specifies how the
!    Voronoi cell generators are to be initialized.
!    0, use the F90 RANDOM_NUMBER routine;
!    1, use the Halton sequence.
!
  integer n
  integer ndim
!
  double precision box_max(ndim)
  double precision box_min(ndim)
  double precision cell_generator(ndim,n)
  integer i
  integer ngen
  integer random_generator
  logical reset
!
  reset = .true.

  do i = 1, n
    call region_sampler ( ndim, box_min, box_max, cell_generator(1,i), &
      random_generator, reset, ngen )
  end do

  return
end
subroutine i_to_halton_vector ( seed, base, ndim, r )
!
!*******************************************************************************
!
!! I_TO_HALTON_VECTOR computes an element of a vector Halton sequence.
!
!
!  Reference:
!
!    J H Halton,
!    Numerische Mathematik,
!    Volume 2, pages 84-90.
!
!  Modified:
!
!    26 February 2001
!
!  Parameters:
!
!    Input, integer SEED, the index of the desired element.
!    Only the absolute value of SEED is considered.  SEED = 0 is allowed,
!    and returns R = 0.
!
!    Input, integer BASE(NDIM), the Halton bases, which should be
!    distinct prime numbers.  This routine only checks that each base
!    is greater than 1.
!
!    Input, integer NDIM, the dimension of the sequence.
!
!    Output, double precision R(NDIM), the SEED-th element of the Halton
!    sequence for the given bases.
!
  integer ndim
!
  integer base(ndim)
  double precision base_inv(ndim)
  integer digit(ndim)
  integer i
  double precision r(ndim)
  integer seed
  integer seed2(ndim)
!
  seed2(1:ndim) = abs ( seed )

  r(1:ndim) = 0.0D+00

  if ( any ( base(1:ndim) <= 1 ) ) then
    write ( *, * ) ' '
    write ( *, * ) 'I_TO_HALTON_VECTOR - Fatal error!'
    write ( *, * ) '  An input base BASE is <= 1!'
    do i = 1, ndim
      write ( *, * ) i, base(i)
    end do
    stop
  end if

  base_inv(1:ndim) = 1.0D+00 / dble ( base(1:ndim) )

  do while ( any ( seed2(1:ndim) /= 0 ) )
    digit(1:ndim) = mod ( seed2(1:ndim), base(1:ndim) )
    r(1:ndim) = r(1:ndim) + dble ( digit(1:ndim) ) * base_inv(1:ndim)
    base_inv(1:ndim) = base_inv(1:ndim) / dble ( base(1:ndim) )
    seed2(1:ndim) = seed2(1:ndim) / base(1:ndim)
  end do

  return
end
subroutine quality ( ndim, n, cell_moment, cell_volume, region_volume )
!
!*******************************************************************************
!
!! QUALITY computes some quality measures for a set of points in a region.
!
!
!  Discussion:
!
!    The quality measures report on how evenly spread the points are.
!
!  Modified:
!
!    02 May 2001
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, integer N, the number of cell generators.
!
!    Input, double precision CELL_MOMENT(NDIM,NDIM,N), the second moment
!    matrix for each Voronoi cell.
!
!    Input, double precision CELL_VOLUME(N), the Voronoi cell volumes.
!
!    Input, double precision REGION_VOLUME, the volume of the region,
!    as input by the user or estimated by VCM.
!
  integer n
  integer ndim
!
  double precision cell_det(n)
  double precision cell_moment(ndim,ndim,n)
  double precision cell_trace(n)
  double precision cell_volume(n)
  double precision dd_l1
  double precision dd_l2
  double precision dd_linf
  double precision dmat_det_2d
  double precision dmat_det_3d
  double precision ev
  double precision ev_l1
  double precision ev_l2
  double precision ev_linf
  integer i
  integer j
  double precision matrix2(2,2)
  double precision matrix3(3,3)
  double precision region_volume
  double precision tr
  double precision tr_l1
  double precision tr_l2
  double precision tr_linf
!
!  Measure 1: the deviation of the cell volumes from the expected cell volume.
!
  ev = region_volume / dble ( n )
  ev_linf = maxval ( abs ( ev - cell_volume(1:n) ) )
  ev_l1 = sum ( abs ( ev - cell_volume(1:n) ) )
  ev_l2 = sqrt ( sum ( ( ev - cell_volume(1:n) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) 'QUALITY'
  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #1:'
  write ( *, * ) '    ( Cell_Volume - Expected Cell Volume )'
  write ( *, * ) ' '
  write ( *, * ) '  Expected Cell Volume = ', ev
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', ev_l1
  write ( *, * ) '  L2 norm =              ', ev_l2
  write ( *, * ) '  L1 norm / N =          ', ev_l1 / dble ( n )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', ev_l2 / sqrt ( dble ( n ) )
  write ( *, * ) '  L-Inf norm =           ', ev_linf
!
!  Measure 2: the deviation of the traces of the cell second moment matrices
!  from the average.
!
  do i = 1, n
    cell_trace(i) = 0.0D+00
    do j = 1, ndim
      cell_trace(i) = cell_trace(i) + cell_moment(j,j,i)
    end do
  end do

  tr = sum ( cell_trace(1:n) ) / dble ( n )
  tr_linf = maxval ( abs ( tr - cell_trace(1:n) ) )
  tr_l1 = sum ( abs ( tr - cell_trace(1:n) ) )
  tr_l2 = sqrt ( sum ( ( tr - cell_trace(1:n) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #2:'
  write ( *, * ) '    ( Cell_Trace - Average Cell Trace )'
  write ( *, * ) ' '
  write ( *, * ) '  Average Cell Trace = ', tr
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', tr_l1
  write ( *, * ) '  L2 norm =              ', tr_l2
  write ( *, * ) '  L1 norm / N =          ', tr_l1 / dble ( n )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', tr_l2 / sqrt ( dble ( n ) )
  write ( *, * ) '  L-Inf norm =           ', tr_linf
!
!  Measure 3: the determinant of the deviatoric matrix
!
  if ( ndim == 2 ) then

    do i = 1, n
      matrix2(1:ndim,1:ndim) = cell_moment(1:ndim,1:ndim,i)
      do j = 1, ndim
        matrix2(j,j) = matrix2(j,j) - cell_trace(i) / dble ( ndim )
      end do

      cell_det(i) = dmat_det_2d ( matrix2 )

    end do

  else if ( ndim == 3 ) then

    do i = 1, n
      matrix3(1:ndim,1:ndim) = cell_moment(1:ndim,1:ndim,i)
      do j = 1, ndim
        matrix3(j,j) = matrix3(j,j) - cell_trace(i) / dble ( ndim )
      end do

      cell_det(i) = dmat_det_3d ( matrix3 )

    end do

  end if

  dd_linf = maxval ( abs ( cell_det(1:n) ) )
  dd_l1 = sum ( abs ( cell_det(1:n) ) )
  dd_l2 = sqrt ( sum ( ( cell_det(1:n) )**2 ) )

  write ( *, * ) ' '
  write ( *, * ) ' '
  write ( *, * ) '  Measure #3:'
  write ( *, * ) '    ( The determinant of the deviatoric matrix )'
  write ( *, * ) ' '
  write ( *, * ) '  L1 norm =              ', dd_l1
  write ( *, * ) '  L2 norm =              ', dd_l2
  write ( *, * ) '  L1 norm / N =          ', dd_l1 / dble ( n )
  write ( *, * ) '  L2 norm / sqrt ( N ) = ', dd_l2 / sqrt ( dble ( n ) )
  write ( *, * ) '  L-Inf norm =           ', dd_linf

  return
end
subroutine region_sampler ( ndim, box_min, box_max, x, random_function, &
  reset, ngen )
!
!******************************************************************************
!
!! REGION_SAMPLER returns a sample point in the physical region.
!
!
!  Discussion:
!
!    The calculations are done in NDIM 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 NDIM is no greater than 3.
!
!  Modified:
!
!    10 April 2001
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision BOX_MIN(NDIM), BOX_MAX(NDIM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, double precision DR, a tolerance used by DIATOM when testing
!    whether a point is within, outside of, or on the boundary of the
!    physical region.
!
!    Output, double precision X(NDIM), the random point.
!
!    Input, integer RANDOM_FUNCTION, specifies the random function.
!    0, uniform random numbers from F90 RANDOM_NUMBER.
!    1, Halton sequence.
!
!    Input/output, logical RESET.
!    If TRUE, the Halton sequence should be reset.
!
!    Input, logical USE_DIATOM, is TRUE if DIATOM is to be called 
!    to determine whether a point lies in the physical region; if it is
!    FALSE than a much simplified routine is used.
!
!    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 ndim
!
  double precision box_max(ndim)
  double precision box_min(ndim)
  double precision density
  integer, save, dimension ( 3 ) :: halton_base = (/ 2, 3, 5 /)
  integer, save :: halton_seed = 1
  integer i
  integer ival
  double precision mdens
  integer ngen
  double precision r(ndim)
  double precision d2
  integer random_function
  logical reset
  double precision x(ndim)
  double precision, parameter :: zero = 0.0D+00
!
  ngen = 0

  if ( reset ) then
    halton_seed = 1
    reset = .false.
  end if

  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 ( *, * ) ' '
      if ( random_function == 0 ) then
        write ( *, * ) '  Using F90 RANDOM_NUMBER.'
      else if ( random_function == 1 ) then
        write ( *, * ) '  Using Halton sequence.'
        write ( *, * ) '  Current seed is ', halton_seed
      end if
      write ( *, * ) ' '
      write ( *, * ) '  Current random value is:'
      write ( *, '(3g14.6)' ) r(1:ndim)
      write ( *, * ) ' '
      write ( *, * ) '  Current sample point is:'
      write ( *, '(3g14.6)' ) x(1:ndim)
      stop
    end if
!
!  Generate a point at random using:
!  0: a uniformly distributed random value;
!  1: a Halton random value.
!
    if ( random_function == 0 ) then

      call random_number ( r(1:ndim) )

    else if ( random_function == 1 ) then

      call i_to_halton_vector ( halton_seed, halton_base, ndim, r )
      halton_seed = halton_seed + 1

    else

      write ( *, * ) ' '
      write ( *, * ) 'REGION_SAMPLER - Fatal error!'
      write ( *, * ) '  Illegal value of RANDOM_FUNCTION = ', random_function
      stop

    end if
!
!  Determine a point in the bounding box.
!
    x(1:ndim) = ( ( 1.0D+00 - r(1:ndim) ) * box_min(1:ndim) &
                            + r(1:ndim)   * box_max(1:ndim) )
!
!  Now determine if the point is in the region.
!
    call test_region ( x, ndim, ival )

    if ( ival == 1 ) then
      exit
    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
!
!  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 test_region ( x, ndim, ival )
!
!*******************************************************************************
!
!! TEST_REGION determines if a point is within the physical region.
!
!
!  Discussion:
!
!    Using a simple routine like this is only appropriate for a simple 
!    region that can be easily defined by user formulas.  This version of 
!    the routine is a demonstration that implements the 2D and 3D versions 
!    of the "tuning fork" test region.
!
!    Computation of the "on-the-boundary" case is not considered important.
!    Only "inside" or "outside" is essential.
!
!  Modified:
!
!    16 April 2001
!
!  Parameters:
!
!    Input, double precision X(NDIM), the point to be checked.
!
!    Input, integer NDIM, the dimension of the space.
!
!    Output, integer IVAL, indicates the status of the point:
!    -1: the point is on the boundary of the region.
!     0: the point is outside the region.
!    +1: the point is inside the region.
!
  integer ndim
!
  double precision c
  integer ival
  double precision x(ndim)
!
  ival = 0

  if ( ndim == 2 ) then

    if ( x(2) >= 40.0D+00 ) then

      if ( x(1) < 45.0D+00 .or. x(1) > 55.0D+00 ) then
        return
      end if

    else

      c = sqrt ( ( x(1) - 50.0D+00 )**2 + x(2)**2 )

      if ( c < 30.0D+00 .or. c > 40.0D+00 ) then
        return
      end if

    end if

  else if ( ndim == 3 ) then

    if ( x(3) <= 20.0D+00 ) then

      if ( x(2) >= 40.0D+00 ) then

        if ( x(1) < 45.0D+00 .or. x(1) > 55.0D+00 ) then
          return
        end if

      else

        c = sqrt ( ( x(1) - 50.0D+00 )**2 + x(2)**2 )

        if ( c < 30.0D+00 .or. c > 40.0D+00 ) then
          return
        end if

      end if

    else

      return

    end if

  end if

  ival = 1

  return
end
subroutine vcm ( ndim, box_min, box_max, n, cell_generator, ns_mom, &
  region_volume_given, region_volume, cell_volume, cell_centroid, cell_moment )
!
!*******************************************************************************
!
!! VCM calculates Voronoi cell volumes, centroids and second moments.
!
!
!  Discussion:
!
!    A Monte Carlo sampling is used to estimate the quantities.
!
!  Modified:
!
!    11 April 2001
!
!  Parameters:
!
!    Input, integer NDIM, the spatial dimension.
!
!    Input, double precision BOX_MIN(NDIM), BOX_MAX(NDIM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, integer N, the number of cell generators.
!
!    Input, double precision CELL_GENERATOR(NDIM,N), the cell generators.
!
!    Input, integer NS_MOM, the number of sample points per cell generator.
!
!    Input, logical REGION_VOLUME_GIVEN,
!    TRUE: the region volume is input in REGION_VOLUME.
!    FALSE: the region volume must be estimated by this routine.
!
!    Input/output, double precision REGION_VOLUME, the volume of the region.
!    If REGION_VOLUME_GIVEN is TRUE, then REGION_VOLUME is input by the user.
!    Otherwise, the volume is estimated and output by this routine.
!
!    Output, double precision CELL_VOLUME(N), the Voronoi cell volumes.
!
!    Output, double precision CELL_CENTROID(NDIM,N), the Voronoi cell centroids.
!
!    Output, double precision CELL_MOMENT(NDIM,NDIM,N), the second moment
!    matrix for each Voronoi cell.
!
  integer ndim
  integer n
!
  double precision box_max(ndim)
  double precision box_min(ndim)
  double precision box_volume
  double precision cell_centroid(ndim,n)
  double precision cell_generator(ndim,n)
  integer cell_hit(n)
  double precision cell_moment(ndim,ndim,n)
  double precision cell_volume(n)
  integer i
  integer j
  integer k
  integer nearest
  integer ngen
  integer ns_mom
  integer ntries
  integer, parameter :: random_mom = 0
  double precision region_volume
  double precision region_volume_estimate
  logical region_volume_given
  logical reset
  double precision x(ndim)
!
!  Zero out the arrays.
!
  cell_centroid(1:ndim,1:n) = 0.0D+00
  cell_hit(1:n) = 0
  cell_moment(1:ndim,1:ndim,1:n) = 0.0D+00
!
!  Sample the region N * NS_MOM times, and keep track of which cell generator
!  is closest to each sampling point.
!
  ntries = 0
  reset = .true.

  do k = 1, n * ns_mom

    call region_sampler ( ndim, box_min, box_max, x, random_mom, reset, ngen )

    ntries = ntries + ngen

    call find_closest ( ndim, x, n, cell_generator, nearest )

    cell_hit(nearest) = cell_hit(nearest) + 1

    cell_centroid(1:ndim,nearest) = cell_centroid(1:ndim,nearest) + x(1:ndim)

    do i = 1, ndim
      do j = 1, ndim
        cell_moment(i,j,nearest) = cell_moment(i,j,nearest) + x(i) * x(j)
      end do
    end do

  end do
!
!  Estimate the area or volume if it was not given.
!
  box_volume = product ( ( box_max(1:ndim)-box_min(1:ndim) ) )
  region_volume_estimate = dble ( n * ns_mom ) * box_volume / dble ( ntries )

  write ( *, * ) ' '
  write ( *, * ) '  Volume of bounding box is     ', box_volume
  if ( region_volume_given ) then
    write ( *, * ) '  Given volume of region is     ', region_volume
  else
    region_volume = region_volume_estimate
  end if
  write ( *, * ) '  Estimated volume of region is ', region_volume_estimate
!
!  Estimate the geometric integrals for each Voronoi cell.
!
  do k = 1, n

    if ( cell_hit(k) > 0 ) then

      cell_volume(k) = dble ( cell_hit(k) ) * region_volume / dble ( n * ns_mom )
      cell_centroid(1:ndim,k) = cell_centroid(1:ndim,k) / dble ( cell_hit(k) )
      cell_moment(1:ndim,1:ndim,k) = cell_moment(1:ndim,1:ndim,k) &
        / dble ( cell_hit(k) )

      do i = 1, ndim
        do j = 1, ndim
          cell_moment(i,j,k) = cell_moment(i,j,k) &
            - cell_centroid(i,k) * cell_centroid(j,k)
        end do
      end do

    end if

  end do

  return
end