subroutine cell_volume_computation ( dim_num, box_min, box_max, cell_num, &
  cell_gen_coord, 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_GEN_COORD(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_gen_coord(dim_num,cell_num)
  real cell_volume(cell_num)
  integer i
  integer j
  integer nearest(sample_num)
  integer ngen
  real region_volume
  real sample_coord(dim_num,sample_num)
!
!  Generate the sampling points.
!
  call region_sampler ( dim_num, sample_num, box_min, box_max, &
    sample_coord, ngen )
!
!  Find closest generators.
!
  call find_closest ( dim_num, sample_num, sample_coord, cell_num, &
    cell_gen_coord, nearest )

  cell_count(1:cell_num) = 0

  do i = 1, sample_num
    j = nearest(i)
    if ( j < 1 .or. j > cell_num ) then
      write ( *, * ) ' '
      write ( *, * ) 'CELL_VOLUME_COMPUTATION - Fatal error!'
      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_dense_iteration ( dim_num, box_min, box_max, cell_num, &
  cell_gen_coord, density_value, density_coord, sample_num )
!
!******************************************************************************
!
!! CVT_DENSE_ITERATION takes one step of the CVT density 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_GEN_COORD(DIM_NUM,CELL_NUM), the Voronoi
!    cell generators.  On output, these have been modified
!
!    Input, real DENSITY_VALUE(CELL_NUM), the density values.
!
!    Input, real DENSITY_COORD(DIM_NUM,CELL_NUM), the coordinates of the
!    density definition points.
!
!    Input, integer SAMPLE_NUM, the number of sample points.
!
  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
  integer cell_count(cell_num)
  real cell_gen_coord(dim_num,cell_num)
  logical, parameter :: debug = .false.
  real density_coord(dim_num,cell_num)
  real density_value(cell_num)
  integer i
  integer nearest(sample_num)
  real sample_coord(dim_num,sample_num)
!
!  Generate sampling points using the current density function.
!
  call region_sampler_weighted ( dim_num, sample_num, cell_num, &
    box_min, box_max, sample_coord, density_value, density_coord, &
    cell_gen_coord )
!
!  For each sample point, find the nearest cell generator.
!
  call find_closest ( dim_num, sample_num, sample_coord, cell_num, &
    cell_gen_coord, nearest )
!
!  Add each sampling point to the averaging data for its nearest generator.
!
  cell_count(1:cell_num) = 0
  cell_gen_coord(1:dim_num,1:cell_num) = 0.0E+00

  write ( *, * ) 'SAMPLE_NUM = ', sample_num
  write ( *, * ) 'CELL_NUM = ', cell_num

  do i = 1, sample_num

    cell = nearest(i)

    if ( cell < 1 .or. cell > cell_num ) then
      write ( *, * ) ' '
      write ( *, * ) '? Illegal value of CELL = ', cell
      write ( *, * ) '  Sample point number is I = ', i
      write ( *, * ) '  Point is ', sample_coord(1:dim_num,i)
    else

      cell_gen_coord(1:dim_num,cell) = cell_gen_coord(1:dim_num,cell) &
        + sample_coord(1:dim_num,i)

      cell_count(cell) = cell_count(cell) + 1

    end if

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

    if ( cell_count(cell) /= 0 ) then

      cell_gen_coord(1:dim_num,cell) = cell_gen_coord(1:dim_num,cell) &
        / real ( cell_count(cell) )

    end if

  end do

  return
end
subroutine find_closest ( dim_num, sample_num, sample_coord, cell_num, &
  cell_gen_coord, nearest )
!
!******************************************************************************
!
!! FIND_CLOSEST finds the Voronoi cell generator closest to a point.
!
!
!  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, integer SAMPLE_NUM, the number of sample points.
!
!    Input, real SAMPLE_COORD(DIM_NUM,SAMPLE_NUM), the points to be checked.
!
!    Input, integer CELL_NUM, the number of cell generatorrs.
!
!    Input, real CELL_GEN_COORD(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
!
  integer cell
  real cell_gen_coord(dim_num,cell_num)
  real distance
  real dist_sq
  integer nearest(sample_num)
  integer sample
  real sample_coord(dim_num,sample_num)
!
  do sample = 1, sample_num

    nearest(sample) = 0
    distance = huge ( distance )

    do cell = 1, cell_num

      dist_sq = sum ( &
        ( cell_gen_coord(1:dim_num,cell) - sample_coord(1:dim_num,sample) )**2 )

      if ( dist_sq < distance ) then
        distance = dist_sq
        nearest(sample) = cell
      end if

    end do

    if ( nearest(sample) == 0 ) then
      write ( *, * ) ' '
      write ( *, * ) 'FIND_CLOSEST - Fatal error!'
      write ( *, * ) '  Sample = ', sample
      write ( *, * ) '  SAMPLE_COORD= ', sample_coord(1:dim_num,sample)
      stop
    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.
!
  implicit none
!
  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_coord, ngen )
!
!******************************************************************************
!
!! REGION_SAMPLER returns sample points 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_COORD(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.
!
  implicit none
!
  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_coord(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_coord(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_coord(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_coord(1:dim_num,j), dim_num, ival )

      if ( ival == 1 ) then
        exit
      end if

    end do

  end do

  return
end
subroutine region_sampler_weighted ( dim_num, sample_num, cell_num, &
  box_min, box_max, sample_coord, density_value, density_coord, &
  cell_gen_coord )
!
!******************************************************************************
!
!! REGION_SAMPLER_WEIGHTED returns weighted sample points
!
!
!  Modified:
!
!    26 September 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer DIM_NUM, the spatial dimension.
!
!    Input, integer SAMPLE_NUM, the number of sample points to compute.
!
!    Input, integer CELL_NUM, the number of Voronoi cells.
!
!    Input, real BOX_MIN(DIM_NUM), BOX_MAX(DIM_NUM), the coordinates
!    of the two extreme corners of the bounding box.
!
!    Input, real DENSITY_VALUE(CELL_NUM), the density values.
!
!    Input, real DENSITY_COORD(DIM_NUM,CELL_NUM), the coordinates of the
!    density definition points.
!
!    Input, real CELL_GEN_COORD(DIM_NUM,CELL_NUM), the coordinates of the
!    Voronoi cell generator points.
!
  implicit none
!
  integer cell_num
  integer dim_num
  integer sample_num
!
  real box_max(dim_num)
  real box_min(dim_num)
  integer cell
  integer cell_count(cell_num)
  real cell_gen_coord(dim_num,cell_num)
  real density_coord(dim_num,cell_num)
  real density_value(cell_num)
  integer i
  integer nearest(1)
  integer ngen
  real r
  integer reject_num
  integer sample
  real sample_coord(dim_num,sample_num)
  real sample_dense
  real sample_dense_prime
!
  reject_num = 0
!
!  Generate sample points according to the current density.
!
  do sample = 1, sample_num

    do
!
!  Pick a uniformly random point in the region.
!
      call region_sampler ( dim_num, 1, box_min, box_max, &
        sample_coord(1,sample), ngen )
!
!  Compute its density.
!
      call spline_linear_val ( cell_num, density_coord, density_value, &
        sample_coord(1,sample), sample_dense, sample_dense_prime )
!
!  Accept the point only if a random number is less than or equal to
!  the density (which is normalized, so that the maximum value is 1).
!
      call random_number ( harvest = r )

      if ( r <= sample_dense ) then
        call find_closest ( dim_num, 1, sample_coord(1,sample), cell_num, &
          cell_gen_coord, nearest )
        cell = nearest(1)
        exit
      else
        reject_num = reject_num + 1
      end if

    end do

  end do

  write ( *, * ) ' '
  write ( *, * ) 'REGION_SAMPLER_WEIGHTED:'
  write ( *, * ) '  Reject = ', reject_num
  write ( *, * ) '  Accept = ', sample_num

  return
end
subroutine test_region ( x, dim_num, 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.
!
!    Computation of the "on-the-boundary" case is not considered important.
!    Only "inside" or "outside" is essential.
!
!  Modified:
!
!    16 April 2001
!
!  Parameters:
!
!    Input, real X(DIM_NUM), the point to be checked.
!
!    Input, integer DIM_NUM, 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.
!
  implicit none
!
  integer dim_num
!
  integer ival
  real x(dim_num)
!
  ival = 0

  if ( dim_num == 1 ) then

    if ( 0.0 <= x(1) .and. x(1) <= 10.0 ) then
      ival = 1
    end if

  else if ( dim_num == 2 ) then

    if ( 0.0 <= x(1) .and. x(1) <= 10.0 .and. &
         0.0 <= x(2) .and. x(2) <= 10.0 ) then
      ival = 1
    end if

  end if

  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 sgtsl ( n, c, d, e, b, info )
!
!*******************************************************************************
!
!! SGTSL solves a general tridiagonal linear system.
!
!
!  Reference:
!
!    Dongarra, Moler, Bunch and Stewart,
!    LINPACK User's Guide,
!    SIAM, (Society for Industrial and Applied Mathematics),
!    3600 University City Science Center,
!    Philadelphia, PA, 19104-2688.
!    ISBN 0-89871-172-X
!
!  Modified:
!
!    31 October 2001
!
!  Parameters:
!
!    Input, integer N, the order of the tridiagonal matrix.
!
!    Input/output, real C(N), contains the subdiagonal of the tridiagonal
!    matrix in entries C(2:N).  On output, C is destroyed.
!
!    Input/output, real D(N).  On input, the diagonal of the matrix.
!    On output, D is destroyed.
!
!    Input/output, real E(N), contains the superdiagonal of the tridiagonal
!    matrix in entries E(1:N-1).  On output E is destroyed.
!
!    Input/output, real B(N).  On input, the right hand side.  On output,
!    the solution.
!
!    Output, integer INFO, error flag.
!    0, normal value.
!    K, the K-th element of the diagonal becomes exactly zero.  The
!       subroutine returns if this error condition is detected.
!
  implicit none
!
  integer n
!
  real b(n)
  real c(n)
  real d(n)
  real e(n)
  integer info
  integer k
  real t
!
  info = 0
  c(1) = d(1)

  if ( n >= 2 ) then

    d(1) = e(1)
    e(1) = 0.0E+00
    e(n) = 0.0E+00

    do k = 1, n - 1
!
!  Find the larger of the two rows, and interchange if necessary.
!
      if ( abs ( c(k+1) ) >= abs ( c(k) ) ) then
        call r_swap ( c(k), c(k+1) )
        call r_swap ( d(k), d(k+1) )
        call r_swap ( e(k), e(k+1) )
        call r_swap ( b(k), b(k+1) )
      end if
!
!  Fail if no nonzero pivot could be found.
!
      if ( c(k) == 0.0E+00 ) then
        info = k
        return
      end if
!
!  Zero elements.
!
      t = -c(k+1) / c(k)
      c(k+1) = d(k+1) + t * d(k)
      d(k+1) = e(k+1) + t * e(k)
      e(k+1) = 0.0E+00
      b(k+1) = b(k+1) + t * b(k)

    end do

  end if

  if ( c(n) == 0.0E+00 ) then
    info = n
    return
  end if
!
!  Back solve.
!
  b(n) = b(n) / c(n)

  if ( n > 1 ) then

    b(n-1) = ( b(n-1) - d(n-1) * b(n) ) / c(n-1)

    do k = n-2, 1, -1
      b(k) = ( b(k) - d(k) * b(k+1) - e(k) * b(k+2) ) / c(k)
    end do

  end if

  return
end
subroutine spline_linear_intset ( int_n, int_x, int_v, data_n, data_x, data_y )
!
!*******************************************************************************
!
!! SPLINE_LINEAR_INTSET sets a linear spline with given integral properties.
!
!
!  Discussion:
!
!    The user has in mind an interval, divided by INT_N+1 points into
!    INT_N intervals.  A linear spline is to be constructed,
!    with breakpoints at the centers of each interval, and extending
!    continuously to the left of the first and right of the last
!    breakpoints.  The constraint on the linear spline is that it is
!    required that it have a given integral value over each interval.
!
!    A tridiagonal linear system of equations is solved for the
!    values of the spline at the breakpoints.
!
!  Modified:
!
!    02 November 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer INT_N, the number of intervals.
!
!    Input, real INT_X(INT_N+1), the points that define the intervals.
!    Interval I lies between INT_X(I) and INT_X(I+1).
!
!    Input, real INT_V(INT_N), the desired value of the integral of the
!    linear spline over each interval.
!
!    Output, integer DATA_N, the number of data points defining the spline.
!    (This is the same as INT_N).
!
!    Output, real DATA_X(DATA_N), DATA_Y(DATA_N), the values of the independent
!    and dependent variables at the data points.  The values of DATA_X are
!    the interval midpoints.  The values of DATA_Y are determined in such
!    a way that the exact integral of the linear spline over interval I
!    is equal to INT_V(I).
!
  implicit none
!
  integer int_n
!
  real c(int_n)
  real d(int_n)
  integer data_n
  real data_x(int_n)
  real data_y(int_n)
  real e(int_n)
  integer info
  real int_v(int_n)
  real int_x(int_n+1)
!
!  Set up the easy stuff.
!
  data_n = int_n
  data_x(1:data_n) = 0.5E+00 * ( int_x(1:data_n) + int_x(2:data_n+1) )
!
!  Set up C, D, E, the coefficients of the linear system.
!
  c(1) = 0.0E+00
  c(2:data_n-1) = 1.0 &
    - ( 0.5 * ( data_x(2:data_n-1) + int_x(2:data_n-1) ) &
    - data_x(1:data_n-2) ) &
    / ( data_x(2:data_n-1) - data_x(1:data_n-2) )
  c(data_n) = 0.0E+00

  d(1) = int_x(2) - int_x(1)

  d(2:data_n-1) = 1.0 &
    + ( 0.5 * ( data_x(2:data_n-1) + int_x(2:data_n-1) ) &
    - data_x(1:data_n-2) ) &
    / ( data_x(2:data_n-1) - data_x(1:data_n-2) ) &
    - ( 0.5 * ( data_x(2:data_n-1) + int_x(3:data_n) ) - data_x(2:data_n-1) ) &
    / ( data_x(3:data_n) - data_x(2:data_n-1) )

  d(data_n) = int_x(data_n+1) - int_x(data_n)

  e(1) = 0.0E+00

  e(2:data_n-1) = ( 0.5 * ( data_x(2:data_n-1) + int_x(3:data_n) ) &
    - data_x(2:data_n-1) ) / ( data_x(3:data_n) - data_x(2:data_n-1) )

  e(data_n) = 0.0E+00
!
!  Set up DATA_Y, which begins as the right hand side of the linear system.
!
  data_y(1) = int_v(1)
  data_y(2:data_n-1) = 2.0E+00 * int_v(2:data_n-1) &
    / ( int_x(3:int_n) - int_x(2:int_n-1) )
  data_y(data_n) = int_v(data_n)
!
!  Solve the linear system.
!
  call sgtsl ( data_n, c, d, e, data_y, info )

  if ( info /= 0 ) then
    write ( *, * ) ' '
    write ( *, * ) 'SPLINE_LINEAR_INTSET - Fatal error!'
    write ( *, * ) '  The linear system is singular.'
    stop
  end if

  return
end
subroutine r_swap ( x, y )
!
!*******************************************************************************
!
!! R_SWAP swaps two real values.
!
!
!  Modified:
!
!    01 May 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, real X, Y.  On output, the values of X and
!    Y have been interchanged.
!
  implicit none
!
  real x
  real y
  real z
!
  z = x
  x = y
  y = z

  return
end
subroutine spline_linear_val ( ndata, tdata, ydata, tval, yval, ypval )
!
!*******************************************************************************
!
!! SPLINE_LINEAR_VAL evaluates a linear spline at a specific point.
!
!
!  Discussion:
!
!    Because of the extremely simple form of the linear spline,
!    the raw data points ( TDATA(I), YDATA(I)) can be used directly to
!    evaluate the spline at any point.  No processing of the data
!    is required.
!
!  Modified:
!
!    06 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer NDATA, the number of data points defining the spline.
!
!    Input, real TDATA(NDATA), YDATA(NDATA), the values of the independent
!    and dependent variables at the data points.  The values of TDATA should
!    be distinct and increasing.
!
!    Input, real TVAL, the point at which the spline is to be evaluated.
!
!    Output, real YVAL, YPVAL, the value of the spline and its first
!    derivative dYdT at TVAL.  YPVAL is not reliable if TVAL is exactly
!    equal to TDATA(I) for some I.
!
  implicit none
!
  integer ndata
!
  integer left
  integer right
  real tdata(ndata)
  real tval
  real ydata(ndata)
  real ypval
  real yval
!
!  Find the interval [ TDATA(LEFT), TDATA(RIGHT) ] that contains, or is
!  nearest to, TVAL.
!
  call rvec_bracket ( ndata, tdata, tval, left, right )
!
!  Now evaluate the piecewise linear function.
!
  ypval = ( ydata(right) - ydata(left) ) / ( tdata(right) - tdata(left) )

  yval = ydata(left) +  ( tval - tdata(left) ) * ypval

  return
end
subroutine rvec_bracket ( n, x, xval, left, right )
!
!*******************************************************************************
!
!! RVEC_BRACKET searches a sorted array for successive brackets of a value.
!
!
!  Discussion:
!
!    If the values in the vector are thought of as defining intervals
!    on the real line, then this routine searches for the interval
!    nearest to or containing the given value.
!
!  Modified:
!
!    06 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, length of input array.
!
!    Input, real X(N), an array sorted into ascending order.
!
!    Input, real XVAL, a value to be bracketed.
!
!    Output, integer LEFT, RIGHT, the results of the search.
!    Either:
!      XVAL < X(1), when LEFT = 1, RIGHT = 2;
!      XVAL > X(N), when LEFT = N-1, RIGHT = N;
!    or
!      X(LEFT) <= XVAL <= X(RIGHT).
!
  implicit none
!
  integer n
!
  integer i
  integer left
  integer right
  real x(n)
  real xval
!
  do i = 2, n - 1

    if ( xval < x(i) ) then
      left = i - 1
      right = i
      return
    end if

   end do

  left = n - 1
  right = n

  return
end
subroutine rvec_sort_heap_a ( n, a )
!
!*******************************************************************************
!
!! RVEC_SORT_HEAP_A ascending sorts a real array using heap sort.
!
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Modified:
!
!    15 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the number of entries in the array.
!
!    Input/output, real A(N).
!    On input, an unsorted array.
!    On output, the array has been sorted.
!
  implicit none
!
  integer n
!
  real a(n)
  integer n1
!
  if ( n <= 1 ) then
    return
  end if
!
!  1: Put A into heap form.
!
  call rvec_heap_a ( n, a )
!
!  2: Sort A.
!
!  The largest object in the heap is in A(1).
!  Move it to position A(N).
!
  call r_swap ( a(1), a(n) )
!
!  Consider the diminished heap of size N1.
!
  do n1 = n-1, 2, -1
!
!  Restore the heap structure of A(1) through A(N1).
!
    call rvec_heap_a ( n1, a )
!
!  Take the largest object from A(1) and move it to A(N1).
!
    call r_swap ( a(1), a(n1) )

  end do

  return
end
subroutine rvec_heap_a ( n, a )
!
!*******************************************************************************
!
!! RVEC_HEAP_A reorders an array of reals into an ascending heap.
!
!
!  Definition:
!
!    An ascending heap is an array A with the property that, for every index J,
!    A(J) >= A(2*J) and A(J) >= A(2*J+1), (as long as the indices
!    2*J and 2*J+1 are legal).
!
!  Diagram:
!
!                  A(1)
!                /      \
!            A(2)         A(3)
!          /     \        /  \
!      A(4)       A(5)  A(6) A(7)
!      /  \       /   \
!    A(8) A(9) A(10) A(11)
!
!  Reference:
!
!    A Nijenhuis and H Wilf,
!    Combinatorial Algorithms,
!    Academic Press, 1978, second edition,
!    ISBN 0-12-519260-6.
!
!  Modified:
!
!    15 April 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer N, the size of the input array.
!
!    Input/output, real A(N).
!    On input, an unsorted array.
!    On output, the array has been reordered into a heap.
!
  implicit none
!
  integer n
!
  real a(n)
  integer i
  integer ifree
  real key
  integer m
!
!  Only nodes N/2 down to 1 can be "parent" nodes.
!
  do i = n/2, 1, -1
!
!  Copy the value out of the parent node.
!  Position IFREE is now "open".
!
    key = a(i)
    ifree = i

    do
!
!  Positions 2*IFREE and 2*IFREE + 1 are the descendants of position
!  IFREE.  (One or both may not exist because they exceed N.)
!
      m = 2 * ifree
!
!  Does the first position exist?
!
      if ( m > n ) then
        exit
      end if
!
!  Does the second position exist?
!
      if ( m + 1 <= n ) then
!
!  If both positions exist, take the larger of the two values,
!  and update M if necessary.
!
        if ( a(m+1) > a(m) ) then
          m = m + 1
        end if

      end if
!
!  If the large descendant is larger than KEY, move it up,
!  and update IFREE, the location of the free position, and
!  consider the descendants of THIS position.
!
      if ( a(m) <= key ) then
        exit
      end if

      a(ifree) = a(m)
      ifree = m

    end do
!
!  Once there is no more shifting to do, KEY moves into the free spot IFREE.
!
    a(ifree) = key

  end do

  return
end