Program e02jdfe ! E02JDF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: e02jdf, nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: liopts = 100, lopts = 100, & nin = 5, nout = 6 ! .. Local Scalars .. Integer :: ifail, lcoefs, lsmaxp, lsminp, & n, nxcels, nycels ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: coefs(:), f(:), x(:), y(:) Real (Kind=nag_wp) :: opts(lopts), pmax(2), pmin(2) Integer :: iopts(liopts) ! .. Intrinsic Procedures .. Intrinsic :: maxval, minval ! .. Executable Statements .. Write (nout,*) 'E02JDF Example Program Results' ! Generate the data to fit. Call generate_data(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs,coefs) ! Initialize the options arrays and set/get some options. Call handle_options(iopts,liopts,opts,lopts) ! Compute the spline coefficients. ifail = 0 Call e02jdf(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs,coefs,iopts,opts, & ifail) ! pmin and pmax form the bounding box of the spline. We must not attempt to ! evaluate the spline outside this box. pmin(:) = (/minval(x),minval(y)/) pmax(:) = (/maxval(x),maxval(y)/) Deallocate (x,y,f) ! Evaluate the approximation at a vector of values. Call evaluate_at_vector(coefs,iopts,opts,pmin,pmax) ! Evaluate the approximation on a mesh. Call evaluate_on_mesh(coefs,iopts,opts,pmin,pmax) Contains Subroutine generate_data(n,x,y,f,lsminp,lsmaxp,nxcels,nycels,lcoefs, & coefs) ! Reads n from a data file and then generates an x and a y vector of n ! pseudorandom uniformly distributed values on (0,1]. These are passed ! to the bivariate function of R. Franke to create the data set to fit. ! The remaining input data for E02JDF are set to suitable values for ! this problem, as discussed by Davydov and Zeilfelder. ! .. Use Statements .. Use nag_library, Only: g05kff, g05saf ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: lseed = 1, mstate = 17 ! .. Scalar Arguments .. Integer, Intent (Out) :: lcoefs, lsmaxp, lsminp, n, & nxcels, nycels ! .. Array Arguments .. Real (Kind=nag_wp), Allocatable, Intent (Out) :: coefs(:), f(:), x(:), & y(:) ! .. Local Scalars .. Integer :: ifail, lstate, subid ! .. Local Arrays .. Integer :: seed(lseed), state(mstate) ! .. Intrinsic Procedures .. Intrinsic :: exp ! .. Executable Statements .. Continue ! Read the size of the data set to be generated and fitted. ! (Skip the heading in the data file.) Read (nin,*) Read (nin,*) n Allocate (x(n),y(n),f(n)) ! Initialize the random number generator and then generate the data. seed(:) = (/32958/) lstate = mstate ifail = 0 Call g05kff(1,subid,seed,lseed,state,lstate,ifail) ifail = 0 Call g05saf(n,state,x,ifail) ifail = 0 Call g05saf(n,state,y,ifail) ! Ensure that the bounding box stretches all the way to (0,0) and (1,1) x(1) = 0.0_nag_wp y(1) = 0.0_nag_wp x(n) = 1.0_nag_wp y(n) = 1.0_nag_wp f(:) = 0.75_nag_wp*exp(-((9._nag_wp*x(:)-2._nag_wp)**2+(9._nag_wp*y(:) & -2._nag_wp)**2)/4._nag_wp) + 0.75_nag_wp*exp(-(9._nag_wp*x(:)+ & 1._nag_wp)**2/49._nag_wp-(9._nag_wp*y(:)+1._nag_wp)/10._nag_wp) + & 0.5_nag_wp*exp(-((9._nag_wp*x(:)-7._nag_wp)**2+(9._nag_wp*y(:)- & 3._nag_wp)**2)/4._nag_wp) - 0.2_nag_wp*exp(-(9._nag_wp*x(:)- & 4._nag_wp)**2-(9._nag_wp*y(:)-7._nag_wp)**2) ! Grid size for the approximation. nxcels = 6 nycels = nxcels ! Identify the computation. Write (nout,*) Write (nout,*) 'Computing the coefficients of a C^1 spline & &approximation to Franke''s function' Write (nout,99999) 'Using a ', nxcels, ' by ', nycels, ' grid' ! Local-approximation control parameters. lsminp = 3 lsmaxp = 100 ! Set up the array to hold the computed spline coefficients. lcoefs = (((nxcels+2)*(nycels+2)+1)/2)*10 + 1 Allocate (coefs(lcoefs)) Return 99999 Format (1X,A,I2,A,I2,A) End Subroutine generate_data Subroutine handle_options(iopts,liopts,opts,lopts) ! Auxiliary routine for initializing the options arrays and ! for demonstrating how to set and get optional parameters. ! .. Use Statements .. Use nag_library, Only: e02zkf, e02zlf ! .. Scalar Arguments .. Integer, Intent (In) :: liopts, lopts ! .. Array Arguments .. Real (Kind=nag_wp), Intent (Out) :: opts(lopts) Integer, Intent (Out) :: iopts(liopts) ! .. Local Scalars .. Real (Kind=nag_wp) :: rvalue Integer :: ifail, ivalue, optype Character (3) :: cvalue Character (80) :: optstr ! .. Executable Statements .. ifail = 0 Call e02zkf('Initialize = E02JDF',iopts,liopts,opts,lopts,ifail) ! Set some non-default parameters for the local approximation method. Write (optstr,99999) 'Minimum Singular Value LPA = ', & 1._nag_wp/32._nag_wp ifail = 0 Call e02zkf(optstr,iopts,liopts,opts,lopts,ifail) ifail = 0 Call e02zkf('Polynomial Starting Degree = 3',iopts,liopts,opts,lopts, & ifail) ! Set some non-default parameters for the global approximation method. ifail = 0 Call e02zkf('Averaged Spline = Yes',iopts,liopts,opts,lopts,ifail) ! As an example of how to get the value of an optional parameter, ! display whether averaging of local approximations is in operation. ifail = 0 Call e02zlf('Averaged Spline',ivalue,rvalue,cvalue,optype,iopts,opts, & ifail) If (cvalue=='YES') Then Write (nout,*) 'Using an averaged local approximation' End If Return 99999 Format (A,E16.9) End Subroutine handle_options Subroutine evaluate_at_vector(coefs,iopts,opts,pmin,pmax) ! Evaluates the approximation at a vector of values. ! .. Use Statements .. Use nag_library, Only: e02jef ! .. Array Arguments .. Real (Kind=nag_wp), Intent (In) :: coefs(*), opts(*), pmax(2), & pmin(2) Integer, Intent (In) :: iopts(*) ! .. Local Scalars .. Integer :: i, ifail, nevalv ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: fevalv(:), xevalv(:), yevalv(:) ! .. Intrinsic Procedures .. Intrinsic :: max, min ! .. Executable Statements .. Read (nin,*) nevalv Allocate (xevalv(nevalv),yevalv(nevalv),fevalv(nevalv)) Read (nin,*)(xevalv(i),yevalv(i),i=1,nevalv) ! Force the points to be within the bounding box of the spline. Do i = 1, nevalv xevalv(i) = max(xevalv(i),pmin(1)) xevalv(i) = min(xevalv(i),pmax(1)) yevalv(i) = max(yevalv(i),pmin(2)) yevalv(i) = min(yevalv(i),pmax(2)) End Do ifail = 0 Call e02jef(nevalv,xevalv,yevalv,coefs,fevalv,iopts,opts,ifail) Write (nout,*) Write (nout,*) 'Values of computed spline at (x_i,y_i):' Write (nout,*) Write (nout,99999) 'x_i', 'y_i', 'f(x_i,y_i)' Write (nout,99998)(xevalv(i),yevalv(i),fevalv(i),i=1,nevalv) Return 99999 Format (1X,3A12) 99998 Format (1X,3F12.2) End Subroutine evaluate_at_vector Subroutine evaluate_on_mesh(coefs,iopts,opts,pmin,pmax) ! Evaluates the approximation on a mesh of n_x * n_y values. ! .. Use Statements .. Use nag_library, Only: e02jff ! .. Implicit None Statement .. Implicit None ! .. Array Arguments .. Real (Kind=nag_wp), Intent (In) :: coefs(*), opts(*), pmax(2), & pmin(2) Integer, Intent (In) :: iopts(*) ! .. Local Scalars .. Integer :: i, ifail, j, nxeval, nyeval Logical :: print_mesh ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: fevalm(:,:), xevalm(:), & yevalm(:) Real (Kind=nag_wp) :: h(2), ll_corner(2), ur_corner(2) ! .. Intrinsic Procedures .. Intrinsic :: max, min, real ! .. Executable Statements .. Read (nin,*) nxeval, nyeval Allocate (xevalm(nxeval),yevalm(nyeval),fevalm(nxeval,nyeval)) ! Define the mesh by its lower-left and upper-right corners. Read (nin,*) ll_corner(1:2) Read (nin,*) ur_corner(1:2) ! Set the mesh spacing and the evaluation points. ! Force the points to be within the bounding box of the spline. h(1) = (ur_corner(1)-ll_corner(1))/real(nxeval-1,nag_wp) h(2) = (ur_corner(2)-ll_corner(2))/real(nyeval-1,nag_wp) Do i = 1, nxeval xevalm(i) = ll_corner(1) + real(i-1,nag_wp)*h(1) xevalm(i) = max(xevalm(i),pmin(1)) xevalm(i) = min(xevalm(i),pmax(1)) End Do Do j = 1, nyeval yevalm(j) = ll_corner(2) + real(j-1,nag_wp)*h(2) yevalm(j) = max(yevalm(j),pmin(2)) yevalm(j) = min(yevalm(j),pmax(2)) End Do ! Evaluate. ifail = 0 Call e02jff(nxeval,nyeval,xevalm,yevalm,coefs,fevalm,iopts,opts,ifail) ! Output the computed function values? Read (nin,*) print_mesh If (.Not. print_mesh) Then Write (nout,*) Write (nout,*) & 'Outputting of the function values on the mesh is disabled' Else Write (nout,*) Write (nout,*) 'Values of computed spline at (x_i,y_j):' Write (nout,*) Write (nout,99999) 'x_i', 'y_j', 'f(x_i,y_j)' Write (nout,99998)((xevalm(i),yevalm(j),fevalm(i, & j),i=1,nxeval),j=1,nyeval) End If Return 99999 Format (1X,3A12) 99998 Format (1X,3F12.2) End Subroutine evaluate_on_mesh End Program e02jdfe