program main c*********************************************************************72 c cc pwl_approx_1d_test() tests pwl_approx_1d(). c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 15 October 2012 c c Author: c c John Burkardt c implicit none integer nc_test_num parameter ( nc_test_num = 4 ) integer nd_test_num parameter ( nd_test_num = 2 ) integer j integer k integer nc integer nc_test(nc_test_num) integer nd integer nd_test(nd_test_num) integer prob integer prob_num save nc_test save nd_test data nc_test / 2, 4, 8, 16 / data nd_test / 16, 64 / call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'pwl_approx_1d_test():' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test pwl_approx_1d().' write ( *, '(a)' ) ' The QR_SOLVE library is needed.' write ( *, '(a)' ) ' The test needs the TEST_INTERP_1D library.' call p00_prob_num ( prob_num ) do prob = 1, prob_num do j = 1, nc_test_num nc = nc_test(j) do k = 1, nd_test_num nd = nd_test(k) call test01 ( prob, nc, nd ) end do end do end do c c Terminate. c write ( *, '(a)' ) '' write ( *, '(a)' ) 'pwl_approx_1d_test():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) return end subroutine test01 ( prob, nc, nd ) c*********************************************************************72 c cc TEST01 tests PWL_APPROX_1D. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 10 October 2012 c c Author: c c John Burkardt c c Parameters: c c Input, integer PROB, the problem index. c c Input, integer NC, the number of control points. c c Input, integer ND, the number of data points. c implicit none integer nc integer nd double precision app_error integer i integer ni integer prob double precision r8vec_norm_affine double precision xc(nc) double precision xd(nd) double precision xi(nd) double precision yc(nc) double precision yd(nd) double precision yi(nd) double precision xmax double precision xmin write ( *, '(a)' ) '' write ( *, '(a)' ) 'TEST01:' write ( *, '(a,i4)' ) & ' Approximate data from TEST_INTERP_1D problem #', prob write ( *, '(a,i4)' ) ' Number of control points = ', nc write ( *, '(a,i4)' ) ' Number of data points = ', nd xmin = 0.0D+00 xmax = 1.0D+00 call r8vec_linspace ( nd, xmin, xmax, xd ) call p00_f ( prob, nd, xd, yd ) if ( nd .lt. 10 ) then call r8vec2_print ( nd, xd, yd, ' Data array:' ) end if c c Determine control values. c call r8vec_linspace ( nc, xmin, xmax, xc ) call pwl_approx_1d ( nd, xd, yd, nc, xc, yc ) c c #1: Does approximant come close to function at data points? c ni = nd do i = 1, ni xi(i) = xd(i) end do call pwl_interp_1d ( nc, xc, yc, ni, xi, yi ) app_error = r8vec_norm_affine ( ni, yi, yd ) / dble ( ni ) write ( *, '(a)' ) '' write ( *, '(a,g14.6)' ) & ' L2 approximation error averaged per data node = ', app_error return end subroutine r8vec_linspace ( n, a, b, x ) c*********************************************************************72 c cc r8vec_linspace() creates a vector of linearly spaced values. c c Discussion: c c An R8VEC is a vector of R8's. c c 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. c c In other words, the interval is divided into N-1 even subintervals, c and the endpoints of intervals are used as the points. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 14 March 2011 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision A, B, the first and last entries. c c Output, double precision X(N), a vector of linearly spaced data. c implicit none integer n double precision a double precision b integer i double precision x(n) if ( n .eq. 1 ) then x(1) = ( a + b ) / 2.0D+00 else do i = 1, n x(i) = ( dble ( n - i ) * a & + dble ( i - 1 ) * b ) & / dble ( n - 1 ) end do end if return end function r8vec_max ( n, a ) c*********************************************************************72 c cc r8vec_max() returns the maximum value in an R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the array. c c Input, double precision A(N), the array. c c Output, double precision R8VEC_MAX, the value of the largest entry. c implicit none integer n double precision a(n) integer i double precision r8_huge parameter ( r8_huge = 1.79769313486231571D+308 ) double precision r8vec_max double precision value value = - r8_huge do i = 1, n value = max ( value, a(i) ) end do r8vec_max = value return end subroutine r8vec2_print ( n, a1, a2, title ) c*********************************************************************72 c cc r8vec2_print() prints an R8VEC2. c c Discussion: c c An R8VEC2 is a dataset consisting of N pairs of R8s, stored c as two separate vectors A1 and A2. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 06 February 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of components of the vector. c c Input, double precision A1(N), A2(N), the vectors to be printed. c c Input, character * ( * ) TITLE, a title. c implicit none integer n double precision a1(n) double precision a2(n) integer i character * ( * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,a,1x,g14.6,2x,g14.6)' ) i, ':', a1(i), a2(i) end do return end subroutine timestamp ( ) c*********************************************************************72 c cc timestamp() prints the YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 June 2014 c c Author: c c John Burkardt c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, & trim ( ampm ) return end