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