program main

!*****************************************************************************80
!
!! shepard_interp_nd_test() tests shepard_interp_nd().
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 June 2024
!
!  Author:
!
!    John Burkardt
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: p_test_num = 4

  integer j
  integer m
  integer n1d
  integer nd
  real ( kind = rk ) p
  real ( kind = rk ), dimension ( p_test_num ) :: p_test = (/ &
    1.0D+00, 2.0D+00, 4.0D+00, 8.0D+00 /)
  integer prob
  integer prob_num

  call timestamp ( )
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'shepard_interp_nd_test():'
  write ( *, '(a)' ) '  Fortran90 version'
  write ( *, '(a)' ) '  Test shepard_interp_nd().'
  write ( *, '(a)' ) '  The R8LIB library is needed.'
  write ( *, '(a)' ) '  This test also needs the TEST_INTERP_ND library.'
!
!  Look at Shepard interpolant on an irregular grid.
!
  nd = 25

  call p00_prob_num ( prob_num )

  do prob = 1, prob_num

    do m = 2, 5, 3

      do j = 1, p_test_num
        p = p_test(j)
        call test01 ( prob, p, m, nd )
      end do

    end do
  end do
!
!  Look at Shepard interpolant on a regular N1D^M grid.
!
  n1d = 5

  call p00_prob_num ( prob_num )

  do prob = 1, prob_num

    do m = 2, 5, 3

      do j = 1, p_test_num
        p = p_test(j)
        call test02 ( prob, p, m, n1d )
      end do

    end do
  end do
!
!  Terminate.
!
  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'shepard_interp_nd_test():'
  write ( *, '(a)' ) '  Normal end of execution.'
  write ( *, '(a)' ) ''
  call timestamp ( )

  stop 0
end
subroutine test01 ( prob, p, m, nd )

!*****************************************************************************80
!
!! TEST01() tests SHEPARD_INTERP() on an irregular grid.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    01 June 2024
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the problem number.
!
!    Input, real ( kind = rk ) P, the power used in the distance weighting.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer ND, the number of data points.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) app_error
  real ( kind = rk ), allocatable :: c(:)
  real ( kind = rk ) int_error
  integer m
  integer nd
  integer ni
  real ( kind = rk ) p
  integer prob
  real ( kind = rk ) r8vec_norm_affine
  real ( kind = rk ), allocatable :: w(:)
  real ( kind = rk ), allocatable :: xd(:,:)
  real ( kind = rk ), allocatable :: xi(:,:)
  real ( kind = rk ), allocatable :: zd(:)
  real ( kind = rk ), allocatable :: ze(:)
  real ( kind = rk ), allocatable :: zi(:)

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'TEST01():'
  write ( *, '(a,i4)' ) '  Interpolate data from TEST_INTERP_ND problem #', prob
  write ( *, '(a,g14.6)' ) '  using Shepard interpolation with P = ', p
  write ( *, '(a,i4)' ) '  spatial dimension M = ', m
  write ( *, '(a,i4,a)' ) '  and an irregular grid of ND = ', nd, ' data points.'
!
!  Set problem parameters:
!
  allocate ( c(1:m) )
  call random_number ( harvest = c(1:m) )
  allocate ( w(1:m) )
  call random_number ( harvest = w(1:m) )

  allocate ( xd(1:m,1:nd) )
  call random_number ( harvest = xd(1:m,1:nd) )

  allocate ( zd(1:nd) )
  call p00_f ( prob, m, c, w, nd, xd, zd )
!
!  #1:  Does interpolant match function at interpolation points?
!
  ni = nd
  allocate ( xi(1:m,1:ni) )
  xi(1:m,1:ni) = xd(1:m,1:ni)
  allocate ( zi(1:ni) )
  call shepard_interp_nd ( m, nd, xd, zd, p, ni, xi, zi )

  int_error = r8vec_norm_affine ( ni, zi, zd ) / real ( ni, kind = rk )

  write ( *, '(a)' ) ''
  write ( *, '(a,g14.6)' ) &
    '  L2 interpolation error averaged per interpolant node = ', int_error

  deallocate ( xi )
  deallocate ( zi )
!
!  #2: Approximation test.  Estimate the integral (f-interp(f))^2.
!
  ni = 1000
  ni = 50
  allocate ( xi(1:m,1:ni) )
  call random_number ( harvest = xi(1:m,1:ni) )

  allocate ( zi(1:ni) )
  call shepard_interp_nd ( m, nd, xd, zd, p, ni, xi, zi )

  allocate ( ze(1:ni) )
  call p00_f ( prob, m, c, w, ni, xi, ze )

  app_error = r8vec_norm_affine ( ni, zi, ze ) / real ( ni, kind = rk )

  write ( *, '(a,g14.6)' ) &
    '  L2 approximation error averaged per 1000 samples =     ', app_error

  deallocate ( c )
  deallocate ( w )
  deallocate ( xd )
  deallocate ( xi )
  deallocate ( zd )
  deallocate ( ze )
  deallocate ( zi )

  return
end
subroutine test02 ( prob, p, m, n1d )

!*****************************************************************************80
!
!! TEST02() tests SHEPARD_INTERP_ND() on a regular N1D^M grid.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    05 August 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer PROB, the problem number.
!
!    Input, real ( kind = rk ) P, the power used in the distance weighting.
!
!    Input, integer M, the spatial dimension.
!
!    Input, integer N1D, the number of points in 1D.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  real ( kind = rk ) a
  real ( kind = rk ) app_error
  real ( kind = rk ) b
  real ( kind = rk ), allocatable :: c(:)
  integer i
  real ( kind = rk ) int_error
  integer m
  integer n1d
  integer nd
  integer ni
  real ( kind = rk ) p
  integer prob
  real ( kind = rk ) r8vec_norm_affine
  real ( kind = rk ), allocatable :: w(:)
  real ( kind = rk ), allocatable :: x1d(:)
  real ( kind = rk ), allocatable :: xd(:,:)
  real ( kind = rk ), allocatable :: xi(:,:)
  real ( kind = rk ), allocatable :: zd(:)
  real ( kind = rk ), allocatable :: ze(:)
  real ( kind = rk ), allocatable :: zi(:)
!
!  Set problem parameters:
!
  allocate ( c(1:m) )
  call random_number ( harvest = c(1:m) )

  allocate ( w(1:m) )
  call random_number ( harvest = w(1:m) )

  nd = n1d ** m

  write ( *, '(a)' ) ''
  write ( *, '(a)' ) 'TEST02():'
  write ( *, '(a,i4)' ) '  Interpolate data from TEST_INTERP_ND problem #', prob
  write ( *, '(a,g14.6)' ) '  using Shepard interpolation with P = ', p
  write ( *, '(a,i4)' ) '  spatial dimension M = ', m
  write ( *, '(a,i6,a)' ) '  and a regular grid of N1D^M = ', nd, ' data points.'

  a = 0.0D+00
  b = 1.0D+00
  allocate ( x1d(1:n1d) )

  call r8vec_linspace ( n1d, a, b, x1d )

  allocate ( xd(1:m,1:nd) )
  do i = 1, m
    call r8vec_direct_product ( i, n1d, x1d, m, nd, xd )
  end do

  allocate ( zd(1:nd) )
  call p00_f ( prob, m, c, w, nd, xd, zd )
!
!  #1:  Does interpolant match function at interpolation points?
!
  ni = nd
  allocate ( xi(1:m,1:nd) )
  xi(1:m,1:nd) = xd(1:m,1:nd)
  allocate ( zi(1:ni) )
  call shepard_interp_nd ( m, nd, xd, zd, p, ni, xi, zi )

  int_error = r8vec_norm_affine ( ni, zi, zd ) / real ( ni, kind = rk )

  write ( *, '(a)' ) ''
  write ( *, '(a,g14.6)' ) &
    '  L2 interpolation error averaged per interpolant node = ', int_error

  deallocate ( xi )
  deallocate ( zi )
!
!  #2: Approximation test.  Estimate the integral (f-interp(f))^2.
!
  ni = 1000
  allocate ( xi(1:m,1:ni) )
  call random_number ( harvest = xi(1:m,1:ni) )

  allocate ( zi(1:ni) )
  call shepard_interp_nd ( m, nd, xd, zd, p, ni, xi, zi )

  allocate ( ze(1:ni) )
  call p00_f ( prob, m, c, w, ni, xi, ze )

  app_error = r8vec_norm_affine ( ni, zi, ze ) / real ( ni, kind = rk )

  write ( *, '(a,g14.6)' ) &
    '  L2 approximation error averaged per 1000 samples =     ', app_error

  deallocate ( c )
  deallocate ( w )
  deallocate ( xd )
  deallocate ( xi )
  deallocate ( zd )
  deallocate ( ze )
  deallocate ( zi )

  return
end
function r8_factorial ( n )

!*****************************************************************************80
!
!! r8_factorial() computes the factorial of N.
!
!  Discussion:
!
!    factorial ( N ) = product ( 1 <= I <= N ) I
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    16 January 1999
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N: the argument of the factorial function.
!    If N is less than 1, the function value is returned as 1.
!
!  Output:
!
!    real ( kind = rk8 ) R8_FACTORIAL: the factorial of N.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) r8_factorial
  integer i
  integer n
  real ( kind = rk8 ) value

  value = 1.0D+00

  do i = 1, n
    value = value * real ( i, kind = rk8 )
  end do

  r8_factorial = value

  return
end
subroutine r8_fake_use ( x )

!*****************************************************************************80
!
!! r8_fake_use() pretends to use an R8 variable.
!
!  Discussion:
!
!    Some compilers will issue a warning if a variable is unused.
!    Sometimes there's a good reason to include a variable in a program,
!    but not to use it.  Calling this function with that variable as
!    the argument will shut the compiler up.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    21 April 2020
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    real ( kind = rk8 ) X, the variable to be "used".
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  real ( kind = rk8 ) x

  if ( x /= x ) then
    write ( *, '(a)' ) '  r8_fake_use(): variable is NAN.'
  end if

  return
end
function r8_mop ( i )

!*****************************************************************************80
!
!! r8_mop() returns the I-th power of -1 as an R8.
!
!  Discussion:
!
!    An R8 is a real ( kind = rk8 ) value.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    07 November 2007
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer I, the power of -1.
!
!  Output:
!
!    real ( kind = rk8 ) R8_MOP, the I-th power of -1.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer i
  real ( kind = rk8 ) r8_mop
  real ( kind = rk8 ) value

  if ( mod ( i, 2 ) == 0 ) then
    value = + 1.0D+00
  else
    value = - 1.0D+00
  end if

  r8_mop = value

  return
end
function r8vec_i4vec_dot_product ( n, r8vec, i4vec )

!*****************************************************************************80
!
!! r8vec_i4vec_dot_product() finds the dot product of an R8VEC and an I4VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    An I4VEC is a vector of I4's.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    30 June 2009
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the dimension of the vectors.
!
!    real ( kind = rk8 ) R8VEC(N), the first vector.
!
!    integer I4VEC(N), the second vector.
!
!  Output:
!
!    real ( kind = rk8 ) R8VEC_I4VEC_DOT_PRODUCT, the dot product.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  integer i4vec(n)
  real ( kind = rk8 ) r8vec(n)
  real ( kind = rk8 ) r8vec_i4vec_dot_product

  r8vec_i4vec_dot_product = dot_product ( r8vec(1:n), &
                                   real ( i4vec(1:n), kind = rk8 ) )

  return
end
subroutine r8vec_direct_product ( factor_index, factor_order, factor_value, &
  factor_num, point_num, x )

!*****************************************************************************80
!
!! r8vec_direct_product() creates a direct product of R8VEC's.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    To explain what is going on here, suppose we had to construct
!    a multidimensional quadrature rule as the product of K rules
!    for 1D quadrature.
!
!    The product rule will be represented as a list of points and weights.
!
!    The J-th item in the product rule will be associated with
!      item J1 of 1D rule 1,
!      item J2 of 1D rule 2,
!      ...,
!      item JK of 1D rule K.
!
!    In particular,
!      X(J) = ( X(1,J1), X(2,J2), ..., X(K,JK))
!    and
!      W(J) = W(1,J1) * W(2,J2) * ... * W(K,JK)
!
!    So we can construct the quadrature rule if we can properly
!    distribute the information in the 1D quadrature rules.
!
!    This routine carries out that task for the abscissas X.
!
!    Another way to do this would be to compute, one by one, the
!    set of all possible indices (J1,J2,...,JK), and then index
!    the appropriate information.  An advantage of the method shown
!    here is that you can process the K-th set of information and
!    then discard it.
!
!  Example:
!
!    Rule 1:
!      Order = 4
!      X(1:4) = ( 1, 2, 3, 4 )
!
!    Rule 2:
!      Order = 3
!      X(1:3) = ( 10, 20, 30 )
!
!    Rule 3:
!      Order = 2
!      X(1:2) = ( 100, 200 )
!
!    Product Rule:
!      Order = 24
!      X(1:24) =
!        ( 1, 10, 100 )
!        ( 2, 10, 100 )
!        ( 3, 10, 100 )
!        ( 4, 10, 100 )
!        ( 1, 20, 100 )
!        ( 2, 20, 100 )
!        ( 3, 20, 100 )
!        ( 4, 20, 100 )
!        ( 1, 30, 100 )
!        ( 2, 30, 100 )
!        ( 3, 30, 100 )
!        ( 4, 30, 100 )
!        ( 1, 10, 200 )
!        ( 2, 10, 200 )
!        ( 3, 10, 200 )
!        ( 4, 10, 200 )
!        ( 1, 20, 200 )
!        ( 2, 20, 200 )
!        ( 3, 20, 200 )
!        ( 4, 20, 200 )
!        ( 1, 30, 200 )
!        ( 2, 30, 200 )
!        ( 3, 30, 200 )
!        ( 4, 30, 200 )
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    18 April 2009
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer FACTOR_INDEX, the index of the factor being
!    processed.  The first factor processed must be factor 1!
!
!    integer FACTOR_ORDER, the order of the factor.
!
!    real ( kind = rk8 ) FACTOR_VALUE(FACTOR_ORDER), the factor values
!    for factor FACTOR_INDEX.
!
!    integer FACTOR_NUM, the number of factors.
!
!    integer POINT_NUM, the number of elements in the
!    direct product.
!
!    real ( kind = rk8 ) X(FACTOR_NUM,POINT_NUM), the elements of
!    the direct product, which are built up gradually.
!
!  Output:
!
!    real ( kind = rk8 ) X(FACTOR_NUM,POINT_NUM), the updated direct product.
!
!  Local:
!
!    integer START, the first location of a block of 
!    values to set.
!
!    integer CONTIG, the number of consecutive values 
!    to set.
!
!    integer SKIP, the distance from the current value 
!    of START to the next location of a block of values to set.
!
!    integer REP, the number of blocks of values to set.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer factor_num
  integer factor_order
  integer point_num

  integer, save :: contig
  integer factor_index
  real ( kind = rk8 ) factor_value(factor_order)
  integer j
  integer k
  integer, save :: rep
  integer, save :: skip
  integer start
  real ( kind = rk8 ) x(factor_num,point_num)

  if ( factor_index == 1 ) then
    contig = 1
    skip = 1
    rep = point_num
    x(1:factor_num,1:point_num) = 0.0D+00
  end if

  rep = rep / factor_order
  skip = skip * factor_order

  do j = 1, factor_order

    start = 1 + ( j - 1 ) * contig

    do k = 1, rep
      x(factor_index,start:start+contig-1) = factor_value(j)
      start = start + skip
    end do

  end do

  contig = contig * factor_order

  return
end
subroutine r8vec_linspace ( n, a, b, x )

!*****************************************************************************80
!
!! r8vec_linspace() creates a vector of linearly spaced values.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
!
!    In other words, the interval is divided into N-1 even subintervals,
!    and the endpoints of intervals are used as the points.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    14 March 2011
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the number of entries in the vector.
!
!    real ( kind = rk8 ) A, B, the first and last entries.
!
!  Output:
!
!    real ( kind = rk8 ) X(N), a vector of linearly spaced data.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) a
  real ( kind = rk8 ) b
  integer i
  real ( kind = rk8 ) x(n)

  if ( n == 1 ) then

    x(1) = ( a + b ) / 2.0D+00

  else

    do i = 1, n
      x(i) = ( real ( n - i,     kind = rk8 ) * a   &
             + real (     i - 1, kind = rk8 ) * b ) &
             / real ( n     - 1, kind = rk8 )
    end do

  end if

  return
end
function r8vec_norm_affine ( n, v0, v1 )

!*****************************************************************************80
!
!! r8vec_norm_affine() returns the affine norm of an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8's.
!
!    The affine vector L2 norm is defined as:
!
!      R8VEC_NORM_AFFINE(V0,V1)
!        = sqrt ( sum ( 1 <= I <= N ) ( V1(I) - V0(I) )^2
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 October 2010
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer N, the order of the vectors.
!
!    real ( kind = rk8 ) V0(N), the base vector.
!
!    real ( kind = rk8 ) V1(N), the vector whose affine norm is desired.
!
!  Output:
!
!    real ( kind = rk8 ) R8VEC_NORM_AFFINE, the L2 norm of V1-V0.
!
  implicit none

  integer, parameter :: rk8 = kind ( 1.0D+00 )

  integer n

  real ( kind = rk8 ) r8vec_norm_affine
  real ( kind = rk8 ) v0(n)
  real ( kind = rk8 ) v1(n)

  r8vec_norm_affine = sqrt ( sum ( ( v0(1:n) - v1(1:n) ) ** 2 ) )

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! timestamp() prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    15 August 2021
!
!  Author:
!
!    John Burkardt
!
  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
  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 ( *, '(i2.2,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 )

  return
end